Registration of Tungsten fibres on XCT images

This demo aims to demonstrate the use of gVirtualXRay and mathematical optimisation to register polygon meshes onto X-ray microtomography (micro-CT) scans of a tungsten fibre. Below is an example of CT slice.

Our simulations include beam-hardening due to polychromatism and they take into account the response of the detector.

We use SimpleGVXR's Python wrapper and Python packages commonly used in tomography reconstruction (Tomopy), image processing (scikit-image and SimpleITK), computer vision (OpenCV), and non-linear numerical optimization (CMA-ES, Covariance Matrix Adaptation Evolution Strategy).

Registration steps

1. Initialisation
   • Import Python packages
   • Global variables with values corresponding to known parameters
   • Load the image data from the experiment at ESRF
   • Recontruct the corresponding CT data
   • Normalise the image data
   • Set the X-ray simulation environment
   • LSF
   • Find circles to identify the centre of fibres
2. Registration of a cube
3. Optimisation of the cores and fibres radii
4. Recentre each core/fibre
5. Optimisation the radii after recentring
6. Optimisation of the beam spectrum
7. Optimisation of the Poisson noise
8. Optimisation of the phase contrast and the radii
9. Optimisation of the phase contrast and the LSF

Import packages

We need to import a few libraries (called packages in Python). We use:

• copy: duplicating images using deepcopies;
• math: the floor function;
• os: creating a new directory
• glob: retrieving file names in a directory;
• numpy: who doesn't use numpy?
• imageio: creating GIF files;
• skimage: comparing the reference CT slice and the simulated one, computing the Radon transform of an image, and perform a CT reconstruction using FBP and SART;
• sklearn: comparing the reference CT slice and the simulated onecomparing the reference CT slice and the simulated one
• tomopy: another package for CT reconstruction;
• SimpleITK: image processing and saving volume data;
• OpenCV (cv2): Hough transform and bilateral filter (an edge-preserving smoothing filter);
• matplotlib: plotting data;
• cma: non-linear numerical optimization;
• lsf: the line spread function to filter the X-ray images; and
• gvxrPython3: simulation of X-ray images using the Beer-Lambert law on GPU.

%matplotlib inline

import copy, math, os, glob

import numpy as np

import imageio
from skimage.transform import iradon
from skimage.util import compare_images
from sklearn.metrics import mean_squared_error
from skimage.util import random_noise

from scipy import ndimage

# import tomopy

import SimpleITK as sitk

import cv2

import cma

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from matplotlib.colors import LogNorm
from matplotlib import cm

plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams['figure.dpi'] = 100 # 200 e.g. is really fine, but slower

from lsf import *

import gvxrPython3 as gvxr

if not os.path.exists("outputs"):
    os.makedirs("outputs");

if not os.path.exists("plots"):
    os.makedirs("plots");

Global variables

We need some global variables.

• g_reference_CT: The reference XCT slice;
• g_reference_sinogram: The Radon transform of the reference XCT slice;
• g_pixel_spacing_in_micrometre and g_pixel_spacing_in_mm: The physical distance between the centre of two successive pixel;
• g_number_of_projections: The total number of angles in the sinogram;
• g_angular_span_in_degrees: The angular span covered by the sinogram;
• g_angular_step: the angular step; and
• g_theta: The rotation angles in degrees (vertical axis of the sinogram).

NoneType = type(None);
pixel_spacing_in_micrometre = 1.9;
pixel_spacing_in_mm = pixel_spacing_in_micrometre * 1e-3;
number_of_projections = 900;
angular_span_in_degrees = 180.0;
angular_step = angular_span_in_degrees / number_of_projections;
theta = np.linspace(0., angular_span_in_degrees, number_of_projections, endpoint=False);
roi_length = 40;

value_range = 6;
num_samples = 15;

bias = None;
gain = None;
scale = None;

fibre_radius = 140 / 2; # um
core_radius = 30 / 2; # um

Load the image data

Load and display the reference projections from a raw binary file, i.e. the target of the registration.

# Target of the registration
reference_normalised_projections = np.fromfile("sino.raw", dtype=np.float32);
reference_normalised_projections.shape = [number_of_projections, int(reference_normalised_projections.shape[0] / number_of_projections)];

def saveMHA(fname, image, spacing):
    """
    save the image into a file.

    :param str fname: the filename
    :param 2D_image image: the image to save
    :param [flt, flt, flt] spacing: the space between two successive voxels along the 3 direction
    """

    #np.savetxt("Tutorial2/outputs/reference_CT.txt", g_reference_CT);
    volume = sitk.GetImageFromArray(image);
    volume.SetSpacing(spacing);
    sitk.WriteImage(volume, fname, useCompression=True);

saveMHA('outputs/reference_normalised_projections.mha', reference_normalised_projections, [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

labels=[theta[0], theta[reference_normalised_projections.shape[0] // 2], theta[-1]];
tics=[0, reference_normalised_projections.shape[0] // 2, reference_normalised_projections.shape[0]-1];
fig=plt.figure();
imgplot = plt.imshow(reference_normalised_projections, cmap="gray");
plt.xlabel("Displacement of projection");
plt.ylabel("Angle of projection (in degrees)");
plt.yticks(tics, labels);
plt.title("Projections after flat-field correction from the experiment at ESRF");
fig.colorbar(imgplot);
plt.savefig('plots/Normalised_projections_from_experiment_ESRF.pdf')
plt.savefig('plots/Normalised_projections_from_experiment_ESRF.png')

In the literature, a projection is often modelled as follows:

$$P = \ln\left(\frac{I_0}{I}\right) = -\ln\left(\frac{I}{I_0}\right) = \sum_n \mu(n) \Delta_x$$

reference_normalised_projections loaded from the binary file corresponds to $\frac{I}{I_0}$. The flat-field correction has already been performed. It is now necessary to linearise the transmission tomography data using:

$$-\ln(normalised\_projections)$$

This new image corresponds to the Radon transform, known as sinogram, of the scanned object in these experimental conditions. Once this is done, we divide the pixels of the sinogram by $\Delta_x$, which is egal to the spacing between two successive pixels along the horizontal axis.

We define a new function to compute the sinogram from flat-field correction and calls it straightaway.

def computeSinogramFromFlatField(normalised_projections):
    """
    This function apply the minus log normalisation
    on the projections that bave been corrected with the flat-field method.

    :param 2D_image normalised_projections: The projections after flat-field corrections
    :return the sinogram.
    """

    # Create a temporary image to hold the sinogram
    simulated_sinogram = copy.deepcopy(normalised_projections);

    # Make sure no value is negative or null (because of the log function)
    # It should not be the case, however, when the Laplacian is used to simulate
    # phase contrast, negative values can be generated.
    threshold = 0.000001
    simulated_sinogram[simulated_sinogram < threshold] = threshold;

    # Apply the minus log normalisation
    simulated_sinogram = -np.log(simulated_sinogram);

    # Rescale the data taking into account the pixel size
    simulated_sinogram /= pixel_spacing_in_micrometre * gvxr.getUnitOfLength("um") / gvxr.getUnitOfLength("cm");

    # Return the new image
    return simulated_sinogram;

reference_sinogram = computeSinogramFromFlatField(reference_normalised_projections);

saveMHA('outputs/reference_sinogram.mha', reference_sinogram, [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

labels=[theta[0], theta[reference_sinogram.shape[0] // 2], theta[-1]];
tics=[0, reference_sinogram.shape[0] // 2, reference_sinogram.shape[0]-1];
fig=plt.figure();
imgplot = plt.imshow(reference_sinogram, cmap="gray");
plt.xlabel("Displacement of projection");
plt.ylabel("Angle of projection (in degrees)");
plt.yticks(tics, labels);
plt.title("Sinogram of the reference image");
fig.colorbar(imgplot);
plt.savefig('plots/Sinogram_reference_image.pdf');
plt.savefig('plots/Sinogram_reference_image.png');

CT reconstruction

Now we got a sinogram, we can reconstruct the CT slice. As we used a synchrotron, we can assume we have a parallel source. It means we can use a FBP rather than the FDK algorithm.

reference_CT = iradon(reference_sinogram.T, theta=theta, circle=True);

saveMHA('outputs/reference_CT.mha', reference_CT, [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

fig=plt.figure();
norm = cm.colors.Normalize(vmax=50, vmin=0)
imgplot = plt.imshow(reference_CT, cmap="gray", norm=norm);
fig.colorbar(imgplot);
plt.title("Reference image (in linear attenuation coefficients, cm$^{-1}$)");
plt.savefig('plots/reference_image_in_mu.pdf');
plt.savefig('plots/reference_image_in_mu.png');

Normalise the image data

Zero-mean unit-variance normalisation is applied to use the reference images in objective functions and perform the registration. Note that it is called standardisation (Z-score Normalisation) in machine learning. It is computed as follows:

$$I' = \frac{I - \bar{I}}{\sigma}$$

Where $I'$ is the image after the original image $I$ has been normalised, $\bar{I}$ is the average pixel value of $I$, and $\sigma$ is its standard deviation.

normalised_reference_sinogram = (reference_sinogram - reference_sinogram.mean()) / reference_sinogram.std();
normalised_reference_CT       = (reference_CT       - reference_CT.mean())       / reference_CT.std();

Set the X-ray simulation environment

First we create an OpenGL context, here using EGL, i.e. no window.

gvxr.createWindow(0, 1, "EGL");
gvxr.setWindowSize(512, 512);

We set the parameters of the X-ray detector (flat pannel), e.g. number of pixels, pixel, spacing, position and orientation:

detector_width_in_pixels = reference_sinogram.shape[1];
detector_height_in_pixels = 1;
distance_object_detector_in_m =    0.08; # = 80 mm

gvxr.setDetectorPosition(-distance_object_detector_in_m, 0.0, 0.0, "m");
gvxr.setDetectorUpVector(0, 1, 0);
gvxr.setDetectorNumberOfPixels(detector_width_in_pixels, detector_height_in_pixels);
gvxr.setDetectorPixelSize(pixel_spacing_in_micrometre, pixel_spacing_in_micrometre, "micrometer");

And the source parameters (beam shape, source position)

# Set up the beam
distance_source_detector_in_m  = 145.0;

gvxr.setSourcePosition(distance_source_detector_in_m - distance_object_detector_in_m,  0.0, 0.0, "m");
gvxr.usePointSource();
# gvxr.useParallelBeam();

The beam spectrum. Here we have a polychromatic beam, with 97% of the photons at 33 keV, 2% at 66 keV and 1% at 99 keV.

energy_spectrum = [(33, 0.97, "keV"), (66, 0.02, "keV"), (99, 0.01, "keV")];

for energy, percentage, unit in energy_spectrum:
    gvxr.addEnergyBinToSpectrum(energy, unit, percentage);

energies_in_keV = [];
weights = [];

for energy, percentage, unit in energy_spectrum:
    weights.append(percentage);
    energies_in_keV.append(energy * gvxr.getUnitOfEnergy(unit) / gvxr.getUnitOfEnergy("keV"));

fig=plt.figure();
plt.xlabel("Energy bin (in keV)");
plt.ylabel("Relative weight");
plt.xticks(energies_in_keV);
plt.yticks(weights);
plt.title("Incident beam spectrum");
plt.bar(energies_in_keV, weights);
plt.savefig('plots/beam_spectrum.pdf');
plt.savefig('plots/beam_spectrum.png');

The material properties (chemical composition and density)

fibre_radius = 140 / 2; # um
fibre_material = [("Si", 0.5), ("C", 0.5)];
fibre_mu = 2.736; # cm-1
fibre_density = 3.2; # g/cm3

core_radius = 30 / 2; # um
core_material = [("W", 1)];
core_mu = 341.61; # cm-1
core_density = 19.3 # g/cm3

g_matrix_width = 0;
g_matrix_height = 0;
g_matrix_x = 0;
g_matrix_y = 0;
matrix_material = [("Ti", 0.9), ("Al", 0.06), ("V", 0.04)];
matrix_mu = 13.1274; # cm-1
matrix_density = 4.42 # g/cm3

The LSF

t = np.arange(-20., 21., 1.);
lsf_kernel=lsf(t*41)/lsf(0);
lsf_kernel/=lsf_kernel.sum();

fig=plt.figure();
plt.title("Response of the detector (LSF)");
plt.plot(t, lsf_kernel);
plt.savefig('plots/LSF.pdf');
plt.savefig('plots/LSF.png');

Find circles to identify the centre of fibres

We can use the Hoguh transform to detect where circles are in the image. However, the input image in OpenCV's function must be in UINT8. We blur it using a bilateral filter (an edge-preserving smoothing filter).

Convert the image to UINT8

We first create a function to convert images in floating point numbers into UINT8.

def float2uint8(anImage):
    uchar_image = copy.deepcopy(anImage);
    uchar_image -= np.min(uchar_image);
    uchar_image /= np.max(uchar_image);
    uchar_image *= 255;
    return uchar_image.astype(np.uint8);

We blur the CT scan using a bilateral filter. It preserves edges.

uint8_reference_CT = float2uint8(reference_CT);
blurred_reference_CT = cv2.bilateralFilter(uint8_reference_CT, 9, 75, 75);

saveMHA('outputs/blurred_reference_CT.mha', blurred_reference_CT, [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

Apply the Hough transform

circles = cv2.HoughCircles(blurred_reference_CT, cv2.HOUGH_GRADIENT, 2, 80,
                            param1=150, param2=5, minRadius=5, maxRadius=15);

Overlay the detected circles on the top of the image

cimg = cv2.cvtColor(blurred_reference_CT, cv2.COLOR_GRAY2BGR);
circles = np.uint16(np.around(circles));

for i in circles[0,:]:

    # draw the outer circle
    cv2.circle(cimg,(i[0],i[1]),i[2],(0,255,0),2);

    # draw the center of the circle
    cv2.circle(cimg,(i[0],i[1]),2,(0,0,255),3);

fig=plt.figure();
imgplot = plt.imshow(cimg);
plt.title("Reference image and detected Tungsten cores");
plt.savefig('plots/fibre_detection_using_Hough_transform.pdf');
plt.savefig('plots/fibre_detection_using_Hough_transform.png');

Unlike the previous example, did did not work that well. Here 13 fibres were missed. Many centres are also misplaced. We will use another technique to register the fibres, the popular Otsu's method. It creates a histogram and uses a heuristic to determine a threshold value.

counts, bins = np.histogram(float2uint8(reference_CT))
plt.hist(bins[:-1], bins, weights=counts)

# Convert the numpy array in float32 into uint8, then into a SITK image
volume = sitk.GetImageFromArray(float2uint8(reference_CT));
volume.SetSpacing([pixel_spacing_in_mm, pixel_spacing_in_mm, pixel_spacing_in_mm]);

# Apply the Otsu's method
otsu_filter = sitk.OtsuThresholdImageFilter();
otsu_filter.SetInsideValue(0);
otsu_filter.SetOutsideValue(1);
seg = otsu_filter.Execute(volume);

# Print the corresponding threshold
print("Threshold:", otsu_filter.GetThreshold());

Threshold: 139.0

sitk.WriteImage(seg, "outputs/cores_segmentation.mha", useCompression=True);

fig=plt.figure();
imgplot = plt.imshow(sitk.GetArrayViewFromImage(sitk.LabelOverlay(volume, seg)));
plt.title("Reference image and detected Tungsten cores");
plt.savefig('plots/fibre_detection_using_otsu_method.pdf');
plt.savefig('plots/fibre_detection_using_otsu_method.png');

Clean up

# Clean-up using mathematical morphology
cleaned_thresh_img = sitk.BinaryOpeningByReconstruction(seg, [3, 3, 3])
cleaned_thresh_img = sitk.BinaryClosingByReconstruction(cleaned_thresh_img, [3, 3, 3])

sitk.WriteImage(cleaned_thresh_img, "outputs/cores_cleaned_segmentation.mha", useCompression=True);

fig=plt.figure();
imgplot = plt.imshow(sitk.GetArrayViewFromImage(sitk.LabelOverlay(volume, cleaned_thresh_img)));
plt.title("Reference image and detected Tungsten cores");
plt.savefig('plots/fibre_detection_using_otsu_method_after_cleaning.pdf');
plt.savefig('plots/fibre_detection_using_otsu_method_after_cleaning.png');

Mark each potential tungsten corewith unique label

core_labels = sitk.ConnectedComponent(cleaned_thresh_img);

fig=plt.figure();
imgplot = plt.imshow(sitk.GetArrayViewFromImage(sitk.LabelOverlay(volume, core_labels)));
plt.title("Cleaned Binary Segmentation of the Tungsten cores");
plt.savefig('plots/fibre_detection_with_label_overlay.pdf');
plt.savefig('plots/fibre_detection_with_label_overlay.png');

Object Analysis

Once we have the segmented objects we look at their shapes and the intensity distributions inside the objects. For each labelled tungsten core, we extract the centroid. Note that sizes and positions are given in millimetres.

shape_stats = sitk.LabelShapeStatisticsImageFilter()
shape_stats.ComputeOrientedBoundingBoxOn()
shape_stats.Execute(core_labels)

centroid_set = [];

for i in shape_stats.GetLabels():
    centroid_set.append(cleaned_thresh_img.TransformPhysicalPointToIndex(shape_stats.GetCentroid(i)));

We now have a list of the centres of all the fibres that can be used as a parameter of the function below to create the cylinders corresponding to the cores and the fibres. For each core, a cylinder is creatd and translated:

gvxr.emptyMesh("core_"  + str(i));
        gvxr.makeCylinder("core_"  + str(i), number_of_sectors, 815.0,  core_radius, "micrometer");
        gvxr.translateNode("core_"  + str(i), y, 0.0, x, "micrometer");

For each fibre, another cylinder is created and translated:

gvxr.emptyMesh("fibre_"  + str(i));
        gvxr.makeCylinder("fibre_"  + str(i), number_of_sectors, 815.0,  fibre_radius, "micrometer");
        gvxr.translateNode("fibre_"  + str(i), y, 0.0, x, "micrometer");

The fibre's cylinder is hollowed to make space for its core:

gvxr.subtractMesh("fibre_" + str(i), "core_" + str(i));

def setFibres(aCentroidSet):
    """This function loads a cylinders in the GPU memory.
    Some are hollow and represent the fibres, some are not and
    correspond to the cores.

    :param array aCentroidSet: a list of cylinder centres.
    """

    global core_radius;
    global fibre_radius;

    # Create empty geometries
    gvxr.emptyMesh("fibre");
    gvxr.emptyMesh("core");

    # Number of sectors to approximate cylinders with triangle meshes
    # It controls the accuracy of the meshes.
    number_of_sectors = 100;

    # Process all the centres from the input list
    for i, cyl in enumerate(aCentroidSet):

        # Convert the centre position from 2D image coordinates in 3D world coordinates
        x = pixel_spacing_in_micrometre * -(cyl[0] - reference_CT.shape[1] / 2 + 0.5);
        y = pixel_spacing_in_micrometre * (cyl[1] - reference_CT.shape[0] / 2 + 0.5);

        # Create empty geometries (is it needed?)
        gvxr.emptyMesh("fibre_" + str(i));
        gvxr.emptyMesh("core_"  + str(i));

        # Create the two corresponding cylinders (fibre and core)
        gvxr.makeCylinder("fibre_" + str(i), number_of_sectors, 815.0, fibre_radius, "micrometer");
        gvxr.makeCylinder("core_"  + str(i), number_of_sectors, 815.0,  core_radius, "micrometer");

        # Translate the two cylinders to the position of their centre
        gvxr.translateNode("fibre_" + str(i), y, 0.0, x, "micrometer");
        gvxr.translateNode("core_"  + str(i), y, 0.0, x, "micrometer");

        # Apply the local transformation matrix (so that we could save the corresponding STL files)
        gvxr.applyCurrentLocalTransformation("fibre_" + str(i));
        gvxr.applyCurrentLocalTransformation("core_" + str(i));

        # Subtract the fibre from the matrix
        gvxr.subtractMesh("matrix", "fibre_" + str(i));

        # Subtract the core from the fibre
        gvxr.subtractMesh("fibre_" + str(i), "core_" + str(i));

        # Save the corresponding STL files
        #gvxr.saveSTLfile("fibre_" + str(i), "Tutorial2/outputs/fibre_" + str(i) + ".stl");
        #gvxr.saveSTLfile("core_" + str(i),  "Tutorial2/outputs/core_"  + str(i) + ".stl");

        # Add the mesh of the current fibre to the overall fibre mesh
        gvxr.addMesh("fibre", "fibre_" + str(i));

        # Add the mesh of the current core to the overall core mesh
        gvxr.addMesh("core",  "core_"  + str(i));

    # Set the mesh colours (for the interactive visualisation)
    gvxr.setColor("fibre", 1.0, 0.0, 0.0, 1.0);
    gvxr.setColor("core",  1.0, 0.0, 1.0, 1.0);

    # Set the fibre's material properties
    #gvxr.setLinearAttenuationCoefficient("fibre", fibre_mu, "cm-1");
    gvxr.setCompound("fibre", "SiC");
    gvxr.setDensity("fibre", fibre_density, "g/cm3");

    # Set the core's material properties
    #gvxr.setLinearAttenuationCoefficient("core", core_mu, "cm-1");
    gvxr.setElement("core", "W");

    # Add the fibres and cores to the X-ray renderer
    gvxr.addPolygonMeshAsInnerSurface("core");
    gvxr.addPolygonMeshAsInnerSurface("fibre");

Registration of a cube

We define a function to create the polygon mesh of the Ti90Al6V4 matrix.

def setMatrix(apGeneSet):
    """This function loads a cube in the GPU memory. The cube represents
    the Ti90Al6V4 matrix.

    apGeneSet[0] is a number between -0.5 and 0.5, related to the translation vector (X component) of the cube. It can be interpreted as a percentage of the detector width.
    apGeneSet[1] is the same as apGeneSet[0], but related to the Y component of the translation vector.
    apGeneSet[2] is a number between -0.5 and 0.5, related to the rotation angle in degrees
    apGeneSet[3] is a scaling factor between -0.5 and 0.5. It can be interpreted as a percentage of the detector width.
    apGeneSet[4] is a scaling factor between -0.5 and 0.5. It can be interpreted as a percentage of apGeneSet[3].
    """

    # Remove all the geometries from the whole scenegraph
    gvxr.removePolygonMeshesFromSceneGraph();

    # Make a cube
    gvxr.makeCube("matrix", 1.0, "micrometer");

    # Translation vector
    x = apGeneSet[0] * detector_width_in_pixels * pixel_spacing_in_micrometre;
    y = apGeneSet[1] * detector_width_in_pixels * pixel_spacing_in_micrometre;
    gvxr.translateNode("matrix", x, 0.0, y, "micrometer");

    # Rotation angle
    rotation_angle_in_degrees = (apGeneSet[2] + 0.5) * 180.0;
    gvxr.rotateNode("matrix", rotation_angle_in_degrees, 0, 1, 0);

    # Scaling factors
    w = (apGeneSet[3] + 0.5) * detector_width_in_pixels * pixel_spacing_in_micrometre;
    h = (apGeneSet[4] + 0.5) * w;
    gvxr.scaleNode("matrix", w, 815, h);

    # Apply the transformation matrix so that we can save the corresponding STL file
    gvxr.applyCurrentLocalTransformation("matrix");

    # Set the matrix's material properties
    gvxr.setMixture("matrix", "Ti90Al6V4");
    gvxr.setDensity("matrix", matrix_density, "g/cm3");

    # Add the matrix to the X-ray renderer
    gvxr.addPolygonMeshAsInnerSurface("matrix");

Simulate the CT acquisition

1. Set the fibre and cores geometries and material properties (Step 39)
2. Set the matrix geometry and material properties (Step 40)
3. Simulate the raw projections for each angle:
   • Without phase contrast (Line 5 of Step 45), or
   • With phase contrast (Lines 14-55 of Step 45)
4. Apply the LSF (Lines 57-60 of Step 45)
5. Apply the flat-field correction (Step 62)
6. Add Poison noise (Step~\ref{??})
7. Apply the minus log normalisation to compute the sinogram (Step 63)

Compute the raw projections and save the data. For this purpose, we define a new function.

def tomographyAcquisition():
    """
    This function simulate a CT acquisition.

    :return the raw projections in keV
    """

    # Crete a new array to save every projection in default unit of energy
    raw_projections = [];

    # For each angle, simulate a projection
    for angle_id in range(0, number_of_projections):

        # Reset the transformation matrix and rotate the scnned object
        gvxr.resetSceneTransformation();
        gvxr.rotateScene(-angular_step * angle_id, 0, 1, 0);

        # Compute the X-ray image
        xray_image = np.array(gvxr.computeXRayImage());

        # Add the projection
        raw_projections.append(xray_image);

    # Convert from the default unit of energy to keV
    raw_projections = np.array(raw_projections);
    raw_projections_in_keV = raw_projections / gvxr.getUnitOfEnergy("keV");

    return raw_projections_in_keV;

Flat-filed correction

Because the data suffers from a fixed-pattern noise in X-ray imaging in actual experiments, it is necessary to perform the flat-field correction of the raw projections using:

$$corrected\_projections = \frac{raw\_projections\_in\_keV − dark\_field\_image}{flat\_field\_image − dark\_field\_image}$$

• $raw\_projections\_in\_keV$ are the raw projections with the X-ray beam turned on and with the scanned object,
• $flat\_field\_image$ is an image with the X-ray beam turned on but without the scanned object, and
• $dark\_field\_image$ is an image with the X-ray beam turned off.

Note that in our example, $raw\_projections\_in\_keV$, $flat\_field\_image$ and $dark\_field\_image$ are in keV whereas $corrected\_projections$ does not have any unit:

$$0 \leq raw\_projections\_in\_keV \leq \sum_E N_0(E) \times E\\0 \leq corrected\_projections \leq 1$$

We define a new function to compute the flat-field correction.

def flatFieldCorrection(raw_projections_in_keV):
    """
    This function applies the flat-field correction on raw projections.

    :param 2D_image raw_projections_in_keV: the raw X-ray projections in keV
    :return the projections (raw_projections_in_keV) after flat-field correction
    """

    # Create a mock dark field image
    dark_field_image = np.zeros(raw_projections_in_keV.shape);

    # Create a mock flat field image
    flat_field_image = np.ones(raw_projections_in_keV.shape);

    # Retrieve the total energy
    total_energy = 0.0;
    energy_bins = gvxr.getEnergyBins("keV");
    photon_count_per_bin = gvxr.getPhotonCountEnergyBins();

    for energy, count in zip(energy_bins, photon_count_per_bin):
        total_energy += energy * count;
    flat_field_image *= total_energy;

    # Apply the actual flat-field correction on the raw projections
    corrected_projections = (raw_projections_in_keV - dark_field_image) / (flat_field_image - dark_field_image);

    return corrected_projections;

The function below is used to simulate a sinogram acquisition. Phase contrast in the projections can be taken into account or not.

def simulateSinogram(sigma_set = None, k_set = None, name_set = None):

    global lsf_kernel;

    # Do not simulate the phase contrast using a Laplacian
    if isinstance(sigma_set, NoneType) or isinstance(k_set, NoneType) or isinstance(name_set, NoneType):

        # Get the raw projections in keV
        raw_projections_in_keV = tomographyAcquisition();

    # Simulate the phase contrast using a Laplacian
    else:

        # Create the convolution filter
        pixel_range = np.linspace(-value_range, value_range, num=int(num_samples), endpoint=True)
        laplacian_kernels = {};

        # Store the L-buffers
        L_buffer_set = {};

        # Look at all the children of the root node
        for label in ["core", "fibre", "matrix"]:
            # Get its L-buffer
            L_buffer_set[label] = getLBuffer(label);

        # Create blank images
        raw_projections_in_keV = np.zeros(L_buffer_set["fibre"].shape);
        phase_contrast_image = np.zeros(L_buffer_set["fibre"].shape);

        for label, k, sigma in zip(name_set, k_set, sigma_set):
            laplacian_kernels[label] = k * laplacian(pixel_range, sigma);

            for z in range(phase_contrast_image.shape[0]):
                for y in range(phase_contrast_image.shape[1]):
                    phase_contrast_image[z][y] += ndimage.convolve((L_buffer_set[label])[z][y], laplacian_kernels[label], mode='wrap');

        for energy, photon_count in zip(gvxr.getEnergyBins("keV"), gvxr.getPhotonCountEnergyBins()):

            # Create a blank image
            attenuation = np.zeros(L_buffer_set["fibre"].shape);

            # Look at all the children of the root node
            #for label in ["core", "fibre", "matrix"]:
            for label in ["core", "fibre", "matrix"]:
                # Get mu for this object for this energy
                mu = gvxr.getLinearAttenuationCoefficient(label, energy, "keV");

                # Compute sum mu * x
                attenuation += L_buffer_set[label] * mu;

            # Store the projection for this energy channel
            raw_projections_in_keV += energy * photon_count * np.exp(-attenuation);

        # Apply the phase contrast
        raw_projections_in_keV -= phase_contrast_image;

    # Apply the LSF line by line
    for z in range(raw_projections_in_keV.shape[0]):
        for y in range(raw_projections_in_keV.shape[1]):
            raw_projections_in_keV[z][y] = ndimage.convolve(raw_projections_in_keV[z][y], lsf_kernel, mode='wrap');

    # Flat-field correction
    normalised_projections = flatFieldCorrection(raw_projections_in_keV);
    normalised_projections[normalised_projections < 0] = 0;

    # Add noise
    if not isinstance(bias, NoneType) and not isinstance(gain, NoneType) and not isinstance(scale, NoneType):

        map = (normalised_projections + (bias + 1)) * gain;
        temp = np.random.poisson(map).astype(np.float);
        temp /= gain;
        temp -= bias + 1;

        # Noise map
        noise_map = (normalised_projections - temp) * scale;
        normalised_projections += noise_map;

    # Linearise
    simulated_sinogram = computeSinogramFromFlatField(normalised_projections);

    return simulated_sinogram, normalised_projections, raw_projections_in_keV;

The function below is used quantify the differences between two images. It is used in the objective function.

def metrics(ref, test):

    normalised_ref = ref.flatten();
    normalised_test = test.flatten();

    normalised_ref  -= normalised_ref.mean();
    normalised_test -= normalised_test.mean();
    normalised_ref  /= normalised_ref.std();
    normalised_test /= normalised_test.std();

    return math.sqrt(mean_squared_error(normalised_ref, normalised_test));

The function below is the objective function used to register the matrix.

def fitnessFunctionCube(x):
    global best_fitness;
    global best_fitness_id;
    global prefix;

    global reference_sinogram;
    global centroid_set;
    global use_fibres;

    global core_radius;
    global fibre_radius;

    # Load the matrix geometrical properties from x
    setMatrix(x);

    # Simulate a sinogram
    simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram();

    # Compute the objective value
    objective = metrics(reference_sinogram, simulated_sinogram);

    # The block below is not necessary for the registration.
    # It is used to save the data to create animations.
    if best_fitness > objective:
        best_fitness = objective;
        matrix_geometry_parameters = copy.deepcopy(x);

        gvxr.saveSTLfile("matrix", "outputs/matrix_" + str(best_fitness_id) + ".stl");

        # Save the simulated sinogram
        simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);
        saveMHA("outputs/" + prefix + "simulated_sinogram_" + str(best_fitness_id) + ".mha",
                simulated_sinogram,
                [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

        # Reconstruct the CT slice
        simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);

        # Save the simulated CT slice
        saveMHA("outputs/" + prefix + "simulated_CT_" + str(best_fitness_id) + ".mha",
                simulated_CT,
                [pixel_spacing_in_mm, pixel_spacing_in_mm, pixel_spacing_in_mm]);

        np.savetxt("outputs/" + prefix + str(best_fitness_id) + ".dat", x, header='x,y,rotation_angle,w,h');

        best_fitness_id += 1;

    return objective

# The registration has already been performed. Load the results.
if os.path.isfile("outputs/cube.dat"):
    matrix_geometry_parameters = np.loadtxt("outputs/cube.dat");
# Perform the registration using CMA-ES
else:
    best_fitness = sys.float_info.max;
    best_fitness_id = 0;
    prefix = "cube_";

    opts = cma.CMAOptions()
    opts.set('tolfun', 1e-2);
    opts['tolx'] = 1e-2;
    opts['bounds'] = [5*[-0.5], 5*[0.5]];

    es = cma.CMAEvolutionStrategy([0.0, 0.0, 0.0, 0.256835938, 0.232903226], 0.5, opts);
    es.optimize(fitnessFunctionCube);

    matrix_geometry_parameters = copy.deepcopy(es.result.xbest);
    np.savetxt("outputs/cube.dat", matrix_geometry_parameters, header='x,y,rotation_angle,w,h');

(4_w,8)-aCMA-ES (mu_w=2.6,w_1=52%) in dimension 5 (seed=769481, Thu Apr 29 17:32:50 2021)
Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]
    1      8 8.797980844271109e-01 1.0e+00 4.38e-01  4e-01  5e-01 0:42.7
    2     16 9.849933135326715e-01 1.2e+00 4.09e-01  3e-01  4e-01 0:47.3
    3     24 9.844200788333163e-01 1.4e+00 3.76e-01  3e-01  4e-01 0:52.2
    4     32 6.930950789572755e-01 1.5e+00 3.61e-01  3e-01  4e-01 1:06.1
    5     40 9.293059250612746e-01 1.7e+00 3.90e-01  3e-01  4e-01 1:10.6
    7     56 9.960254989946028e-01 2.0e+00 3.85e-01  3e-01  4e-01 1:19.7
    9     72 7.597052043204545e-01 2.6e+00 3.76e-01  2e-01  4e-01 1:28.9
   11     88 6.867311644411925e-01 2.4e+00 3.59e-01  2e-01  4e-01 1:47.3
   13    104 6.544822080250594e-01 2.4e+00 2.96e-01  2e-01  3e-01 2:05.3
   15    120 5.707478194537363e-01 2.7e+00 2.72e-01  1e-01  3e-01 2:34.0
   16    128 4.985816567613370e-01 2.6e+00 2.45e-01  1e-01  3e-01 2:47.6
   18    144 4.695495935288412e-01 2.6e+00 1.91e-01  9e-02  2e-01 3:07.0
   21    168 6.178429671801369e-01 3.0e+00 1.81e-01  7e-02  2e-01 3:20.6
   22    176 3.926868606385446e-01 3.2e+00 1.58e-01  6e-02  2e-01 3:43.4
   25    200 5.538982728201552e-01 3.8e+00 1.30e-01  5e-02  1e-01 4:03.4
   29    232 3.946832477630928e-01 3.3e+00 1.12e-01  3e-02  1e-01 4:23.0
   32    256 3.643712812001308e-01 4.1e+00 1.04e-01  3e-02  9e-02 4:50.4
   34    272 3.213786647100584e-01 3.9e+00 8.34e-02  2e-02  7e-02 5:08.5
   37    296 2.692508523848967e-01 4.4e+00 6.27e-02  1e-02  5e-02 5:30.4
   42    336 2.825299071714703e-01 4.1e+00 3.43e-02  7e-03  2e-02 5:52.9
   44    352 2.649162054596793e-01 4.0e+00 3.35e-02  6e-03  2e-02 6:19.3
   46    368 2.555880181704406e-01 4.0e+00 2.78e-02  5e-03  2e-02 6:45.5
   49    392 2.512765648152919e-01 4.4e+00 2.29e-02  3e-03  1e-02 7:07.9
   51    408 2.457158073206211e-01 4.7e+00 1.58e-02  2e-03  9e-03 7:34.0

Apply the result of the registration

# Save the result
setMatrix(matrix_geometry_parameters);
gvxr.saveSTLfile("matrix", "outputs/matrix.stl");

# Translation vector
x = matrix_geometry_parameters[0] * detector_width_in_pixels * pixel_spacing_in_micrometre;
y = matrix_geometry_parameters[1] * detector_width_in_pixels * pixel_spacing_in_micrometre;

# Rotation angle
rotation_angle_in_degrees = (matrix_geometry_parameters[2] + 0.5) * 180.0;

# Scaling factors
w = (matrix_geometry_parameters[3] + 0.5) * detector_width_in_pixels * pixel_spacing_in_micrometre;
h = (matrix_geometry_parameters[4] + 0.5) * w;

print("Matrix");
print("\tposition:", x, y, "um");
print("\trotation:", rotation_angle_in_degrees, "deg");
print("\tsize:", w, h, "um");

Matrix
	position: -55.82766230642614 -22.80317233285621 um
	rotation: 90.78122652385588 deg
	size: 1454.1285834091495 1073.8384786733918 um

Simulate the correspond CT acquisition

# Simulate a sinogram
simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram();
simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);

# Reconstruct the CT slice
simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);
normalised_simulated_CT = (simulated_CT - simulated_CT.mean()) / simulated_CT.std();

# Compute the ZNCC
print("ZNCC matrix registration:",
      "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, normalised_simulated_CT))));

ZNCC matrix registration: 45.29

Display the result of the registration as an animation

def createAnimation(aPrefix, anOutputFile):
    # Find all the images from the output directory
    files = sorted(
        glob.glob(aPrefix + "[0-9]*.mha"))

    # Store the images
    registration_image_set = [];

    # Create the GIF file
    with imageio.get_writer(anOutputFile, mode='I') as writer:

        # Store the PNG filenames
        png_filename_set = [];

        # Process all the images
        for i in range(len(files)):

            # Create the filenames
            mha_fname = aPrefix + str(i) + ".mha";
            png_filename_set.append(aPrefix + str(i) + ".png");

            # Open the MHA file
            float_image = sitk.ReadImage(mha_fname);

            # Convert in a Numpy array
            narray = sitk.GetArrayFromImage(float_image);

            narray -= narray.mean();
            narray /= narray.std();
            registration_image_set.append(narray);

            # Create the figure
            fig, (ax1, ax2, ax3) = plt.subplots(1, 3)

            # Dispay the reference, registration and error map
            fig.suptitle('Registration: Result ' + str(i+1) + "/" + str(len(files)))
            plt.tight_layout();
            norm = cm.colors.Normalize(vmax=1.25, vmin=-0.5)

            # Reference
            ax1.set_title("Reference image");
            imgplot1 = ax1.imshow(normalised_reference_CT, cmap="gray",
                                 norm=norm);

            # Registration
            ax2.set_title("Simulated CT slice after automatic registration");
            imgplot2 = ax2.imshow(narray,
                                 cmap='gray',
                                 norm=norm);

            # Error map
            comp_equalized = compare_images(normalised_reference_CT, narray, method='checkerboard');
            ax3.set_title("Checkboard comparison between\nthe reference and simulated images\nZNCC: " + "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, narray))));
            imgplot3 = ax3.imshow(comp_equalized,
                                 cmap='gray',
                                 norm=norm);

            plt.tight_layout();

            # Save the figure as a PNG file
            plt.savefig(png_filename_set[i])

            # Close the figure
            plt.close()

            # Open the PNG file with imageio and add it to the GIF file
            image = imageio.imread(png_filename_set[i])
            writer.append_data(image)

            # Delete the PNG file
            os.remove(png_filename_set[i]);

        for i in range(15):
            writer.append_data(image)

    return registration_image_set;

if not os.path.exists("plots/cube_registration.gif"):
    cube_registration_image_set = createAnimation("outputs/cube_simulated_CT_",
                'plots/cube_registration.gif');

Adding the fibres

The radius of a tungsten core is 30 / 2 um. The pixel spacing is 1.9 um. The radius in number of pixels is $15/1.9 \approx 7.89$. The area of a core is $(15/1.9)^2 \pi \approx 196$ pixels.

setMatrix(matrix_geometry_parameters);
setFibres(centroid_set);

# Simulate a sinogram
simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram();

simulated_sinogram.shape     = (simulated_sinogram.size     // simulated_sinogram.shape[2],     simulated_sinogram.shape[2]);
normalised_projections.shape = (normalised_projections.size // normalised_projections.shape[2], normalised_projections.shape[2]);
raw_projections_in_keV.shape = (raw_projections_in_keV.size // raw_projections_in_keV.shape[2], raw_projections_in_keV.shape[2]);

saveMHA("outputs/simulated_sinogram_with_fibres.mha",
        simulated_sinogram,
        [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

saveMHA("outputs/normalised_projections_with_fibres.mha",
        normalised_projections,
        [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

saveMHA("outputs/raw_projections_in_keV_with_fibres.mha",
        raw_projections_in_keV,
        [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

# Reconstruct the CT slice
simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);
normalised_simulated_CT = (simulated_CT - simulated_CT.mean()) / simulated_CT.std();

# Compute the ZNCC
print("ZNCC matrix registration with fibres:",
      "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, normalised_simulated_CT))));

ZNCC matrix registration with fibres: 56.03

saveMHA("outputs/simulated_CT_with_fibres.mha",
        simulated_CT,
        [pixel_spacing_in_mm, pixel_spacing_in_mm, pixel_spacing_in_mm]);

saveMHA("outputs/normalised_simulated_CT_with_fibres.mha",
        normalised_simulated_CT,
        [pixel_spacing_in_mm, pixel_spacing_in_mm, pixel_spacing_in_mm]);

norm = cm.colors.Normalize(vmax=1.25, vmin=-0.5)

fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
plt.tight_layout()
fig.suptitle('CT slice with fibres after the registration of the matrix')

ax1.set_title("Reference image");
imgplot1 = ax1.imshow(normalised_reference_CT, cmap="gray",
                     norm=norm);

ax2.set_title("Simulated CT slice after automatic registration");
imgplot2 = ax2.imshow(normalised_simulated_CT,
                     cmap='gray',
                     norm=norm);

comp_equalized = compare_images(normalised_reference_CT, normalised_simulated_CT, method='checkerboard');
ax3.set_title("Checkboard comparison between\n" +
              "the reference and simulated images\nZNCC: " +
              "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, normalised_simulated_CT))));
imgplot3 = ax3.imshow(comp_equalized,
                     cmap='gray',
                     norm=norm);

plt.savefig('plots/simulated_CT_slice_with_fibres_after_cube_registration.pdf');
plt.savefig('plots/simulated_CT_slice_with_fibres_after_cube_registration.png');

Optimisation of the cores and fibres radii

The function below is the objective function used to optimise the radii of the cores and fibres.

def fitnessFunctionFibres(x):
    global best_fitness;
    global best_fitness_id;
    global fibre_radius;
    global core_radius;
    global prefix;

    # Get the radii
    fibre_radius = x[0];
    core_radius = fibre_radius * x[1];

    # Load the matrix
    setMatrix(matrix_geometry_parameters);

    # Load the cores and fibres
    setFibres(centroid_set);

    # Simulate a sinogram
    simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram();

    # Compute the objective value
    objective = metrics(reference_sinogram, simulated_sinogram);

    # The block below is not necessary for the registration.
    # It is used to save the data to create animations.
    if best_fitness > objective:
        best_fitness = objective;

        gvxr.saveSTLfile("core",  "outputs/" + prefix + str(best_fitness_id) + "_cores.stl");
        gvxr.saveSTLfile("fibre", "outputs/" + prefix + str(best_fitness_id) + "_fibres.stl");

        # Save the simulated sinogram
        simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);
        saveMHA("outputs/" + prefix + "simulated_sinogram_" + str(best_fitness_id) + ".mha",
                simulated_sinogram,
                [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

        # Reconstruct the CT slice
        simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);

        # Save the simulated CT slice
        saveMHA("outputs/" + prefix + "simulated_CT_" + str(best_fitness_id) + ".mha",
                simulated_CT,
                [pixel_spacing_in_mm, pixel_spacing_in_mm, pixel_spacing_in_mm]);

        np.savetxt("outputs/" + prefix + str(best_fitness_id) + ".dat", x, header='x,y,rotation_angle,w,h');

        best_fitness_id += 1;

    return objective

# The registration has already been performed. Load the results.
if os.path.isfile("outputs/fibre1_radii.dat"):
    temp = np.loadtxt("outputs/fibre1_radii.dat");
    core_radius = temp[0];
    fibre_radius = temp[1];
# Perform the registration using CMA-ES
else:
    ratio = core_radius / fibre_radius;

    x0 = [fibre_radius, ratio];
    bounds = [[5, 0.01], [1.5 * fibre_radius, 0.95]];

    best_fitness = sys.float_info.max;
    best_fitness_id = 0;
    prefix = "fibre1_";

    opts = cma.CMAOptions()
    opts.set('tolfun', 1e-3);
    opts['tolx'] = 1e-3;
    opts['bounds'] = bounds;

    es = cma.CMAEvolutionStrategy(x0, 0.9, opts);
    es.optimize(fitnessFunctionFibres);
    fibre_radius = es.result.xbest[0];
    core_radius = fibre_radius * es.result.xbest[1];

    np.savetxt("outputs/fibre1_radii.dat", [core_radius, fibre_radius], header='core_radius_in_um,fibre_radius_in_um');

(3_w,6)-aCMA-ES (mu_w=2.0,w_1=63%) in dimension 2 (seed=870684, Thu Apr 29 17:40:56 2021)
Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]
    1      6 3.125542712911615e-01 1.0e+00 8.21e-01  7e-01  8e-01 0:27.3
    2     12 3.791611724817031e-01 1.4e+00 1.30e+00  1e+00  1e+00 0:37.0
    3     18 2.428265395643834e-01 1.4e+00 1.60e+00  1e+00  2e+00 0:55.5
    4     24 3.835705743030111e-01 1.3e+00 1.70e+00  1e+00  2e+00 1:05.2
    5     30 3.917382638265611e-01 1.7e+00 2.01e+00  1e+00  3e+00 1:14.9
    6     36 3.196300202390601e-01 1.8e+00 2.08e+00  2e+00  2e+00 1:24.7
    7     42 2.373196893222082e-01 1.6e+00 1.95e+00  1e+00  2e+00 1:43.0
    8     48 3.421039258375037e-01 1.6e+00 1.78e+00  1e+00  2e+00 1:53.0
    9     54 2.537170288591002e-01 1.7e+00 1.97e+00  1e+00  2e+00 2:02.7
   10     60 2.903177233957678e-01 1.2e+00 1.46e+00  9e-01  1e+00 2:12.4
   12     72 3.447008918873501e-01 1.6e+00 1.54e+00  9e-01  2e+00 2:31.9
   14     84 3.223356324917725e-01 2.0e+00 1.49e+00  8e-01  1e+00 2:51.4
   16     96 2.785877543879386e-01 2.1e+00 1.28e+00  6e-01  9e-01 3:10.8
   18    108 2.985976983941393e-01 1.6e+00 9.75e-01  4e-01  6e-01 3:34.2
   19    114 2.404322674312792e-01 1.5e+00 7.93e-01  3e-01  5e-01 3:48.7
   21    126 3.779575035602915e-01 2.1e+00 7.47e-01  2e-01  5e-01 4:13.0
   23    138 3.785752317751608e-01 2.2e+00 1.40e+00  4e-01  9e-01 4:34.6
   25    150 3.791709273418579e-01 2.3e+00 1.59e+00  5e-01  1e+00 4:54.4
   27    162 2.776782537952786e-01 1.7e+00 9.59e-01  3e-01  5e-01 5:13.9
   29    174 3.776323582191883e-01 2.7e+00 6.51e-01  1e-01  4e-01 5:33.4
   32    192 2.436858734832649e-01 3.3e+00 6.32e-01  9e-02  5e-01 6:02.6
   35    210 2.604841235640152e-01 4.4e+00 3.90e-01  5e-02  2e-01 6:31.8
   38    228 2.430088641059025e-01 3.0e+00 3.33e-01  5e-02  1e-01 7:01.0
   41    246 2.432987885726698e-01 3.6e+00 2.12e-01  2e-02  5e-02 7:30.2
   44    264 2.425078007077185e-01 3.0e+00 1.33e-01  1e-02  2e-02 7:59.4
   47    282 2.422829136657573e-01 2.1e+00 9.70e-02  8e-03  1e-02 8:28.7
   48    288 2.424456091661499e-01 1.7e+00 7.14e-02  5e-03  7e-03 8:38.4

if not os.path.exists("plots/fibre1_registration.gif"):
    registration_image_set = createAnimation("outputs/fibre1_simulated_CT_",
                'plots/fibre1_registration.gif');

Apply the result of the registration

# Load the matrix
setMatrix(matrix_geometry_parameters);

# Load the cores and fibres
setFibres(centroid_set);

gvxr.saveSTLfile("fibre", "outputs/fibre1_fibre.stl");
gvxr.saveSTLfile("core",  "outputs/fibre1_core.stl");

print("Core diameter:", round(core_radius * 2), "um");
print("Fibre diameter:", round(fibre_radius * 2), "um");

# Simulate the corresponding CT aquisition
simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram();
simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);

# Reconstruct the CT slice
simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);
normalised_simulated_CT = (simulated_CT - simulated_CT.mean()) / simulated_CT.std();

# Compute the ZNCC
print("ZNCC radii registration:",
      "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, normalised_simulated_CT))));

Core diameter: 12 um
Fibre diameter: 146 um
ZNCC radii registration: 78.70

The 3D view of the registration looks like:

Recentre each core/fibre

Each fibre is extracted from both the reference CT slice and simulated CT slice. The displacement between the corresponding fibres is computed to maximise the ZNCC between the two. The centre of the fibre is then adjusted accordingly.

def refineCentrePositions(centroid_set, reconstruction_CT_fibres):

    # Exhaustive local search to refine the centre of each cylinder
    roi_length = 40;
    new_centroid_set = [];
    for i, cyl in enumerate(centroid_set):

        centre = [
            cyl[0],
            cyl[1]
        ];

        # extract ROI from reference image
        reference_image = copy.deepcopy(reference_CT[centre[1] - roi_length:centre[1] + roi_length, centre[0] - roi_length:centre[0] + roi_length]);

        # Normalise ROI
        reference_image -= reference_image.mean();
        reference_image /= reference_image.std();

        best_ZNCC = -1;
        best_x_offset = 0;
        best_y_offset = 0;

        for y in range(-10, 11):
            for x in range(-10, 11):

                centre = [
                    cyl[0] + x,
                    cyl[1] + y
                ];

                # extract ROI from test image
                test_image = copy.deepcopy(reconstruction_CT_fibres[centre[1] - roi_length:centre[1] + roi_length, centre[0] - roi_length:centre[0] + roi_length]);

                # Normalise ROI
                test_image -= test_image.mean();
                test_image /= test_image.std();

                # Compare the ROIs
                zncc = np.mean(np.multiply(reference_image.flatten(), test_image.flatten()));

                if best_ZNCC < zncc:
                    best_ZNCC = zncc;
                    best_x_offset = x;
                    best_y_offset = y;

        # Correct the position of the centre of the fibre
        new_centroid_set.append([cyl[0] - best_x_offset, cyl[1] - best_y_offset]);

    return new_centroid_set;

centroid_set = refineCentrePositions(centroid_set, normalised_simulated_CT);

Applying the result of recentring

# Load the matrix
setMatrix(matrix_geometry_parameters);

# Load the cores and fibres
setFibres(centroid_set);

gvxr.saveSTLfile("fibre", "outputs/fibre2_fibre.stl");
gvxr.saveSTLfile("core",  "outputs/fibre2_core.stl");

# Simulate the corresponding CT aquisition
simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram();
simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);

# Reconstruct the CT slice
simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);
normalised_simulated_CT = (simulated_CT - simulated_CT.mean()) / simulated_CT.std();

# Compute the ZNCC
print("ZNCC recentring registration:",
      "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, normalised_simulated_CT))));

ZNCC recentring registration: 82.84

Optimisation the radii after recentring

After recentring the centres, another run of optimisation is executed to refine the radii of the fibres and cores.

# The registration has already been performed. Load the results.
if os.path.isfile("outputs/fibre3_radii.dat"):
    temp = np.loadtxt("outputs/fibre3_radii.dat");
    core_radius = temp[0];
    fibre_radius = temp[1];
# Perform the registration using CMA-ES
else:
    ratio = core_radius / fibre_radius;

    x0 = [fibre_radius, ratio];
    bounds = [[5, 0.01], [1.5 * fibre_radius, 0.95]];

    best_fitness = sys.float_info.max;
    best_fitness_id = 0;
    prefix = "fibre3_";

    opts = cma.CMAOptions()
    opts.set('tolfun', 1e-2);
    opts['tolx'] = 1e-2;
    opts['bounds'] = bounds;

    es = cma.CMAEvolutionStrategy(x0, 0.9, opts);
    es.optimize(fitnessFunctionFibres);
    fibre_radius = es.result.xbest[0];
    core_radius = fibre_radius * es.result.xbest[1];

    np.savetxt("outputs/fibre3_radii.dat", [core_radius, fibre_radius], header='core_radius_in_um,fibre_radius_in_um');

(3_w,6)-aCMA-ES (mu_w=2.0,w_1=63%) in dimension 2 (seed=864832, Thu Apr 29 17:49:59 2021)
Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]
    1      6 2.302303201438057e-01 1.0e+00 9.52e-01  8e-01  1e+00 0:36.5
    2     12 2.246419728555043e-01 1.5e+00 1.41e+00  1e+00  2e+00 0:55.1
    3     18 3.294125411862520e-01 1.6e+00 1.28e+00  1e+00  1e+00 1:04.8
    4     24 2.791682464055642e-01 1.5e+00 1.70e+00  1e+00  3e+00 1:14.6
    5     30 3.468742819162343e-01 2.5e+00 1.60e+00  1e+00  2e+00 1:24.3
    6     36 3.749544429567822e-01 2.3e+00 2.02e+00  2e+00  3e+00 1:34.0
    7     42 3.211507831079822e-01 2.3e+00 2.50e+00  3e+00  4e+00 1:43.8
    8     48 3.465821159554363e-01 2.7e+00 2.73e+00  3e+00  4e+00 1:53.7
    9     54 4.102625895890153e-01 2.4e+00 2.92e+00  3e+00  5e+00 2:03.4
   10     60 3.373677298097274e-01 2.0e+00 3.65e+00  3e+00  7e+00 2:13.2
   12     72 2.756656947845372e-01 2.5e+00 2.37e+00  2e+00  4e+00 2:32.7
   14     84 4.165891843131491e-01 2.3e+00 2.82e+00  2e+00  4e+00 2:52.2
   16     96 3.006621205634746e-01 2.7e+00 2.59e+00  2e+00  3e+00 3:11.7
   18    108 4.100312385036549e-01 4.2e+00 4.73e+00  3e+00  9e+00 3:31.1
   20    120 2.617895428807701e-01 3.5e+00 3.56e+00  2e+00  6e+00 3:50.7
   22    132 3.780213508499477e-01 2.8e+00 3.00e+00  2e+00  4e+00 4:10.2
   24    144 3.734303565927052e-01 3.2e+00 1.90e+00  8e-01  2e+00 4:29.7
   26    156 3.679219432253409e-01 5.3e+00 1.94e+00  7e-01  3e+00 4:49.2
   28    168 3.705810345179621e-01 6.1e+00 1.47e+00  4e-01  2e+00 5:08.7
   30    180 2.994333750485400e-01 3.9e+00 1.69e+00  6e-01  2e+00 5:28.2
   33    198 3.824050556371876e-01 3.1e+00 2.14e+00  8e-01  2e+00 5:57.4
   36    216 2.194526759932435e-01 3.6e+00 2.26e+00  7e-01  2e+00 6:35.3
   38    228 2.161837968266745e-01 4.9e+00 4.08e+00  1e+00  5e+00 7:03.6
   41    246 1.723472792589521e-01 5.0e+00 4.69e+00  1e+00  5e+00 7:41.8
   44    264 1.363911503303434e-01 4.5e+00 8.47e+00  2e+00  1e+01 8:24.2
   46    276 1.218164363740912e-01 4.7e+00 9.96e+00  2e+00  1e+01 8:54.0
   49    294 2.095788515707647e-01 4.3e+00 7.33e+00  1e+00  5e+00 9:23.2
   52    312 1.261216272985952e-01 5.9e+00 8.54e+00  1e+00  7e+00 9:52.4
   55    330 1.272108236926771e-01 7.0e+00 5.30e+00  7e-01  4e+00 10:21.8
   58    348 1.472179769499971e-01 9.4e+00 5.23e+00  5e-01  4e+00 10:56.8
   61    366 1.333231202904199e-01 8.6e+00 2.98e+00  2e-01  2e+00 11:29.2
   65    390 1.392799496086632e-01 8.1e+00 3.22e+00  3e-01  2e+00 12:08.0
   69    414 1.667691906975682e-01 1.1e+01 3.03e+00  1e-01  2e+00 12:47.0
   73    438 1.388294446868860e-01 2.2e+01 1.45e+00  3e-02  7e-01 13:25.9
   77    462 1.322588568119216e-01 7.2e+01 2.39e+00  3e-02  2e+00 14:05.2
   81    486 1.241788531148540e-01 1.4e+02 4.11e+00  3e-02  4e+00 14:44.4
   83    498 1.143051984158120e-01 2.1e+02 4.73e+00  3e-02  6e+00 15:21.8
   86    516 1.172158398153264e-01 4.2e+02 4.10e+00  2e-02  5e+00 16:00.1
   89    534 1.142360988280998e-01 4.1e+02 2.83e+00  1e-02  3e+00 16:47.2
   92    552 1.134408242319352e-01 5.3e+02 1.39e+00  4e-03  1e+00 17:34.2
   95    570 1.133142278427942e-01 5.3e+02 1.02e+00  2e-03  6e-01 18:30.2
   99    594 1.133118257771092e-01 4.7e+02 4.61e-01  7e-04  2e-01 19:26.9
  100    600 1.133133140013076e-01 4.1e+02 3.85e-01  5e-04  1e-01 19:36.9
  101    606 1.133114445843006e-01 3.9e+02 3.38e-01  4e-04  1e-01 19:55.7

if not os.path.exists("plots/fibre3_registration.gif"):
    registration_image_set = createAnimation("outputs/fibre3_simulated_CT_",
                'plots/fibre3_registration.gif');

Apply the result of the registration

# Load the matrix
setMatrix(matrix_geometry_parameters);

# Load the cores and fibres
setFibres(centroid_set);

gvxr.saveSTLfile("fibre", "outputs/fibre3_fibre.stl");
gvxr.saveSTLfile("core",  "outputs/fibre3_core.stl");

print("Core diameter:", round(core_radius * 2), "um");
print("Fibre diameter:", round(fibre_radius * 2), "um");

Core diameter: 15 um
Fibre diameter: 105 um

# Simulate the corresponding CT aquisition
simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram();
simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);

# Reconstruct the CT slice
simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);
normalised_simulated_CT = (simulated_CT - simulated_CT.mean()) / simulated_CT.std();

# Compute the ZNCC
print("ZNCC radii registration:",
      "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, normalised_simulated_CT))));

ZNCC radii registration: 88.82

Optimisation of the Poisson noise

def fitnessFunctionNoise(x):
    global best_fitness;
    global best_fitness_id;
    global prefix;

    bias = x[0];
    gain = x[1];
    scale = x[2];

    # Poisson noise
    map = (normalised_projections_ROI + (bias + 1)) * gain;
    temp = np.random.poisson(map).astype(np.float);
    temp /= gain;
    temp -= bias + 1;

    # Noise map
    noise_map = normalised_projections_ROI - temp;
    noise_map *= scale;
    noisy_image = normalised_projections_ROI + noise_map;

    # Compute the standard deviation of the pixel values in the ROI extracted from the simulated image with noise
    noisy_image_noise_ROI_stddev = 0;
    for y in range(noisy_image.shape[0]):
        noisy_image_noise_ROI_stddev += noisy_image[y].std();
    noisy_image_noise_ROI_stddev /= noisy_image.shape[0];

    # Difference of std dev between the reference and the simulated image
    diff = reference_noise_ROI_stddev - noisy_image_noise_ROI_stddev;
    objective = diff * diff;

    # The block below is not necessary for the registration.
    # It is used to save the data to create animations.
    if best_fitness > objective:
        best_fitness = objective;

        # Save the simulated CT slice
        saveMHA("outputs/" + prefix + "noisy_image_" + str(best_fitness_id) + ".mha",
                noisy_image,
                [pixel_spacing_in_mm, pixel_spacing_in_mm, pixel_spacing_in_mm]);

        np.savetxt("outputs/" + prefix + str(best_fitness_id) + ".dat", [bias, gain, scale], header='bias, gain, scale');

        best_fitness_id += 1;

    return objective

# The registration has already been performed. Load the results.
if os.path.isfile("outputs/poisson-noise.dat"):
    temp = np.loadtxt("outputs/poisson-noise.dat");
    bias = temp[0];
    gain = temp[1];
    scale = temp[2];

# Perform the registration using CMA-ES
else:

    # Extract a ROI from the reference where no object is
    reference_noise_ROI = copy.deepcopy(reference_normalised_projections[450:550,0:125]);

    saveMHA("outputs/reference_noise_ROI.mha",
           reference_noise_ROI,
           [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

    # Compute the standard deviation of the pixel values in the ROI extracted from the reference
    reference_noise_ROI_stddev = 0;
    for y in range(reference_noise_ROI.shape[0]):
        reference_noise_ROI_stddev += reference_noise_ROI[y].std();
    reference_noise_ROI_stddev /= reference_noise_ROI.shape[0];

    # Copy the simulated projection in a temporary variable
    temp = copy.deepcopy(normalised_projections);
    temp.shape = reference_normalised_projections.shape

    # Extract the corresponding ROI
    normalised_projections_ROI = temp[450:550,0:125];

    saveMHA("outputs/normalised_projections_ROI.mha",
           normalised_projections_ROI,
           [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

    # Initialise the values
    bias = 0.0;
    gain = 255.0;
    scale = 1;

    x0 = [bias, gain, scale];

    bounds = [
        [-1.0,   0.0, 0.0],
        [ 5.0, 255.0, 255.0]
    ];

    opts = cma.CMAOptions()
    opts.set('tolfun', 1e-8);
    opts['tolx'] = 1e-8;
    opts['bounds'] = bounds;
    opts['CMA_stds'] = [1, 10, 10];

    best_fitness = sys.float_info.max;
    best_fitness_id = 0;
    prefix = "poisson-noise_";

    es = cma.CMAEvolutionStrategy(x0, 0.25, opts);
    es.optimize(fitnessFunctionNoise);

    bias = es.result.xbest[0];
    gain = es.result.xbest[1];
    scale = es.result.xbest[2];

    np.savetxt("outputs/poisson-noise.dat", [bias, gain, scale], header='bias, gain, scale');

(3_w,7)-aCMA-ES (mu_w=2.3,w_1=58%) in dimension 3 (seed=798878, Thu Apr 29 18:10:17 2021)
Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]
    1      7 2.399548485783231e-04 1.0e+00 1.99e-01  2e-01  2e+00 0:00.0
    2     14 2.397132862885106e-04 1.2e+00 2.10e-01  2e-01  2e+00 0:00.0
    3     21 3.914303339231155e-03 1.5e+00 2.38e-01  3e-01  2e+00 0:00.1
   57    399 3.319192661012088e-11 4.6e+03 3.73e-02  8e-04  8e-01 0:01.0

print("Noise parameters: ", bias, gain, scale)

Noise parameters:  0.7028384906355933 254.93875951097365 0.04895933795291745

Apply the result of the optimisation

# Simulate the corresponding CT aquisition
simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram();
simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);

# Reconstruct the CT slice
simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);
normalised_simulated_CT = (simulated_CT - simulated_CT.mean()) / simulated_CT.std();

# Compute the ZNCC
print("ZNCC noise registration:",
      "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, normalised_simulated_CT))));

ZNCC noise registration: 87.55

norm = cm.colors.Normalize(vmax=1.25, vmin=-0.5)

fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
plt.tight_layout()
fig.suptitle('CT slice with fibres after the registration of the matrix')

ax1.set_title("Reference image");
imgplot1 = ax1.imshow(normalised_reference_CT, cmap="gray",
                     norm=norm);

ax2.set_title("Simulated CT slice after automatic registration");
imgplot2 = ax2.imshow(normalised_simulated_CT,
                     cmap='gray',
                     norm=norm);

comp_equalized = compare_images(normalised_reference_CT, normalised_simulated_CT, method='checkerboard');
ax3.set_title("Checkboard comparison between\n" +
              "the reference and simulated images\nZNCC: " +
              "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, normalised_simulated_CT))));
imgplot3 = ax3.imshow(comp_equalized,
                     cmap='gray',
                     norm=norm);

plt.savefig('plots/simulated_CT_slice_with_noise.pdf');
plt.savefig('plots/simulated_CT_slice_with_noise.png');

Optimisation of the beam spectrum

def fitnessHarmonics(x):

    global energy_spectrum;

    global best_fitness;
    global best_fitness_id;
    global prefix;

    energy_33_keV = x[0];
    first_order_harmonics = x[1];
    second_order_harmonics = x[2];

    # Normalise the beam spectrum
    total = energy_33_keV + first_order_harmonics + second_order_harmonics;
    energy_33_keV /= total;
    first_order_harmonics /= total;
    second_order_harmonics /= total;

    # The beam specturm. Here we have a polychromatic beam.
    gvxr.resetBeamSpectrum();
    energy_spectrum = [(33, energy_33_keV, "keV"), (66, first_order_harmonics, "keV"), (99, second_order_harmonics, "keV")];

    for energy, percentage, unit in energy_spectrum:
        gvxr.addEnergyBinToSpectrum(energy, unit, percentage);

    # Simulate a sinogram
    simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram();

    # Compute the objective value (no normalisation here)
    objective = math.sqrt(mean_squared_error(reference_sinogram.flatten(), simulated_sinogram.flatten()));

    # The block below is not necessary for the registration.
    # It is used to save the data to create animations.
    if best_fitness > objective:
        best_fitness = objective;

        # Save the simulated sinogram
        simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);
        saveMHA("outputs/" + prefix + "simulated_sinogram_" + str(best_fitness_id) + ".mha",
                simulated_sinogram,
                [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

        # Reconstruct the CT slice
        simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);

        # Save the simulated CT slice
        saveMHA("outputs/" + prefix + "simulated_CT_" + str(best_fitness_id) + ".mha",
                simulated_CT,
                [pixel_spacing_in_mm, pixel_spacing_in_mm, pixel_spacing_in_mm]);

        np.savetxt("outputs/" + prefix + str(best_fitness_id) + ".dat", np.array(x) / total, header='33keV,66keV,99keV');

        best_fitness_id += 1;

    return objective

# The registration has already been performed. Load the results.
if os.path.isfile("outputs/spectrum.dat"):
    temp = np.loadtxt("outputs/spectrum.dat");

    # The beam specturm. Here we have a polychromatic beam.
    energy_spectrum = [(33, temp[0], "keV"), (66, temp[1], "keV"), (99, temp[2], "keV")];

# Perform the registration using CMA-ES
else:
    ratio = core_radius / fibre_radius;

    x0 = [0.97, 0.2, 0.1];
    bounds = [[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]];

    best_fitness = sys.float_info.max;
    best_fitness_id = 0;
    prefix = "spectrum_";

    opts = cma.CMAOptions()
    opts.set('tolfun', 1e-4);
    opts['tolx'] = 1e-4;
    opts['bounds'] = bounds;

    es = cma.CMAEvolutionStrategy(x0, 0.25, opts);
    es.optimize(fitnessHarmonics);

    total = es.result.xbest[0] + es.result.xbest[1] + es.result.xbest[2];
    energy_spectrum = [(33, es.result.xbest[0] / total, "keV"), (66, es.result.xbest[1] / total, "keV"), (99, es.result.xbest[2] / total, "keV")];

    np.savetxt("outputs/spectrum.dat", [es.result.xbest[0] / total, es.result.xbest[1] / total, es.result.xbest[2] / total], header='weight of main energy,weight of first order harmonics,weight of second order harmonics');

(3_w,7)-aCMA-ES (mu_w=2.3,w_1=58%) in dimension 3 (seed=851494, Thu Apr 29 18:10:29 2021)
Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]
    1      7 6.790233777334277e+02 1.0e+00 2.36e-01  2e-01  2e-01 0:36.5
    2     14 2.710266191591077e+03 1.4e+00 2.28e-01  2e-01  2e-01 0:46.6
    3     21 1.618025960597630e+03 1.6e+00 2.26e-01  2e-01  2e-01 0:56.7
    4     28 1.716806461023481e+03 1.6e+00 2.52e-01  2e-01  3e-01 1:06.8
    5     35 2.571945696029042e+03 1.9e+00 2.19e-01  1e-01  2e-01 1:16.9
    6     42 5.759469255536430e+02 1.9e+00 1.97e-01  1e-01  2e-01 1:35.8
    7     49 4.486008918813712e+02 1.9e+00 2.01e-01  1e-01  2e-01 1:54.7
    8     56 8.190725918544050e+02 2.2e+00 1.95e-01  9e-02  2e-01 2:04.9
    9     63 5.417669289158210e+02 2.6e+00 1.53e-01  6e-02  2e-01 2:15.0
   10     70 4.174401023648967e+02 2.9e+00 1.32e-01  5e-02  1e-01 2:33.8
   11     77 4.394674174980944e+02 3.0e+00 1.13e-01  4e-02  1e-01 2:43.9
   12     84 3.953122264633439e+02 3.2e+00 1.12e-01  4e-02  1e-01 3:02.8
   14     98 4.999718542793537e+02 3.7e+00 1.97e-01  6e-02  2e-01 3:23.0
   16    112 3.945675298668301e+02 4.4e+00 2.38e-01  7e-02  4e-01 3:52.0
   17    119 3.900960817736810e+02 6.3e+00 2.60e-01  7e-02  4e-01 4:10.8
   19    133 3.928230591120026e+02 6.5e+00 2.31e-01  6e-02  3e-01 4:31.1
   21    147 3.907277089483744e+02 6.6e+00 1.88e-01  5e-02  2e-01 4:51.3
   23    161 3.948360156423612e+02 7.3e+00 1.61e-01  4e-02  2e-01 5:11.5
   25    175 3.966271839726701e+02 9.0e+00 1.52e-01  3e-02  2e-01 5:32.1
   27    189 5.242566304307235e+02 7.6e+00 1.43e-01  2e-02  2e-01 5:52.6
   29    203 3.935636266659243e+02 8.1e+00 1.11e-01  2e-02  1e-01 6:13.1
   31    217 3.897618327805572e+02 9.8e+00 8.84e-02  1e-02  1e-01 6:42.4
   34    238 3.913801325418385e+02 1.2e+01 7.05e-02  7e-03  8e-02 7:13.1
   37    259 3.897304612321485e+02 1.6e+01 3.77e-02  3e-03  4e-02 7:52.6
   39    273 3.894933262406266e+02 2.2e+01 4.48e-02  4e-03  5e-02 8:21.9
   42    294 3.895780376780284e+02 3.1e+01 2.75e-02  3e-03  2e-02 8:52.6
   45    315 3.894829577444672e+02 3.4e+01 2.11e-02  2e-03  2e-02 9:23.4
   47    329 3.894331567245317e+02 3.6e+01 1.86e-02  2e-03  1e-02 9:52.7
   49    343 3.892485509085648e+02 3.8e+01 2.16e-02  2e-03  2e-02 10:30.5
   52    364 3.888720698202041e+02 5.6e+01 3.38e-02  4e-03  3e-02 11:18.3
   55    385 3.889063369477800e+02 6.0e+01 3.39e-02  3e-03  3e-02 11:57.9
   59    413 3.890537671431833e+02 9.0e+01 2.27e-02  2e-03  2e-02 12:38.7
   63    441 3.890194362192212e+02 1.3e+02 1.87e-02  9e-04  1e-02 13:19.5
   67    469 3.889461764107460e+02 1.5e+02 1.01e-02  4e-04  7e-03 14:00.3
   71    497 3.888825909095691e+02 1.8e+02 8.15e-03  2e-04  6e-03 14:41.0
   75    525 3.887492573392220e+02 1.4e+02 1.33e-02  4e-04  8e-03 15:21.8
   79    553 3.889611977751908e+02 1.7e+02 1.44e-02  5e-04  9e-03 16:02.7
   83    581 3.888369837402643e+02 2.0e+02 1.17e-02  3e-04  6e-03 16:43.4
   87    609 3.887892204275699e+02 2.2e+02 1.30e-02  3e-04  7e-03 17:24.4
   91    637 3.888313116581915e+02 1.7e+02 8.74e-03  2e-04  4e-03 18:05.3
   95    665 3.888332698420073e+02 8.8e+01 8.07e-03  1e-04  3e-03 18:46.1
  100    700 3.888865237609597e+02 7.6e+01 9.57e-03  1e-04  3e-03 19:37.1
  105    735 3.888605559460909e+02 8.2e+01 1.07e-02  1e-04  3e-03 20:28.0
  110    770 3.888985450769374e+02 1.1e+02 5.59e-03  6e-05  1e-03 21:19.0
  115    805 3.887237074550530e+02 1.8e+02 4.14e-03  4e-05  1e-03 22:10.0
  119    833 3.889466314584904e+02 2.0e+02 3.10e-03  3e-05  5e-04 22:59.6
  124    868 3.888564845540565e+02 2.3e+02 4.48e-03  4e-05  1e-03 23:50.8
  128    896 3.888922611424571e+02 2.6e+02 4.10e-03  4e-05  1e-03 24:40.5
  133    931 3.889666470088185e+02 2.6e+02 4.39e-03  4e-05  1e-03 25:31.5
  138    966 3.888200746888679e+02 1.7e+02 5.37e-03  4e-05  1e-03 26:22.5
  143   1001 3.888261272611612e+02 2.6e+02 4.34e-03  3e-05  9e-04 27:14.4
  148   1036 3.887948599991024e+02 2.8e+02 2.37e-03  1e-05  4e-04 28:12.8
  152   1064 3.886768883068419e+02 2.4e+02 2.09e-03  1e-05  3e-04 29:23.3
  157   1099 3.888711450168244e+02 2.8e+02 2.27e-03  1e-05  3e-04 30:23.0
  163   1141 3.888457577722545e+02 2.7e+02 2.26e-03  9e-06  2e-04 31:24.4
  169   1183 3.889049897856885e+02 4.5e+02 1.75e-03  7e-06  2e-04 32:25.0
  175   1225 3.888909149612533e+02 3.7e+02 1.11e-03  3e-06  1e-04 33:25.6

if not os.path.exists("plots/spectrum_registration.gif"):
    registration_image_set = createAnimation("outputs/spectrum_simulated_CT_",
                'plots/spectrum_registration.gif');

# Apply the result of the registration
gvxr.resetBeamSpectrum();
for energy, percentage, unit in energy_spectrum:
    gvxr.addEnergyBinToSpectrum(energy, unit, percentage);

for channel in energy_spectrum:
    print(channel);

(33, 0.9565783895953861, 'keV')
(66, 0.01737075206037302, 'keV')
(99, 0.026050858344240937, 'keV')

# Simulate the corresponding CT aquisition
simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram();
simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);

# Reconstruct the CT slice
simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);
normalised_simulated_CT = (simulated_CT - simulated_CT.mean()) / simulated_CT.std();

# Compute the ZNCC
print("ZNCC spectrum registration:",
      "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, normalised_simulated_CT))));

ZNCC spectrum registration: 87.59

Optimisation of the phase contrast and the radii

def laplacian(x, sigma):
    """
    This function create a Laplacian kernel with

    $$ g''(x) = \left(\frac{x^2}{\sigma^4} - \frac{1}{\sigma^2}\right) \exp\left(-\frac{x^2}{2\sigma^2}\right) $$
    
    :param array x: 
    :param float sigma:
    :return the convolution kernel
    """

    return (np.power(x, 2.) / math.pow(sigma, 4) - 1. / math.pow(sigma, 2)) * np.exp(-np.power(x, 2.) / (2. * math.pow(sigma, 2)));

def getLBuffer(object):

    """
    This function compute the L-buffer of the object over all the angles
    
    :param str object: the name of the object 
    :return the L-buffer over all the angles
    """

    # An empty L-buffer
    L_buffer = [];

    # Get the line of L-buffer for each angle
    for angle_id in range(0, number_of_projections):
        gvxr.resetSceneTransformation();
        gvxr.rotateScene(-angular_step * angle_id, 0, 1, 0);

        # Compute the X-ray image
        line_of_L_buffer = np.array(gvxr.computeLBuffer(object));

        # Add the projection
        L_buffer.append(line_of_L_buffer);

    # Return as a numpy array
    return np.array(L_buffer);

def fitnessFunctionLaplacian(x):
    global best_fitness;
    global best_fitness_id;
    global prefix;

    global fibre_radius;
    global core_radius;

    sigma_core = x[0];
    k_core = x[1];
    sigma_fibre = x[2];
    k_fibre = x[3];
    sigma_matrix = x[4];
    k_matrix = x[5];
    core_radius = x[6];
    fibre_radius = x[7];

    # Load the matrix
    setMatrix(matrix_geometry_parameters);

    # Load the cores and fibres
    setFibres(centroid_set);

    # Simulate a sinogram
    simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram(
        [sigma_core, sigma_fibre, sigma_matrix],
        [k_core, k_fibre, k_matrix],
        ["core", "fibre", "matrix"]
    );

    # Compute the objective value
    objective = metrics(reference_sinogram, simulated_sinogram);

    # The block below is not necessary for the registration.
    # It is used to save the data to create animations.
    if best_fitness > objective:
        best_fitness = objective;

        # Save the simulated sinogram
        simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);
        saveMHA("outputs/" + prefix + "simulated_sinogram_" + str(best_fitness_id) + ".mha",
                simulated_sinogram,
                [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

        # Reconstruct the CT slice
        simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);

        # Save the simulated CT slice
        saveMHA("outputs/" + prefix + "simulated_CT_" + str(best_fitness_id) + ".mha",
                simulated_CT,
                [pixel_spacing_in_mm, pixel_spacing_in_mm, pixel_spacing_in_mm]);

        np.savetxt("outputs/" + prefix + str(best_fitness_id) + ".dat", [sigma_core, k_core, sigma_fibre, k_fibre, sigma_matrix, k_matrix, core_radius, fibre_radius], header='sigma_core, k_core, sigma_fibre, k_fibre, sigma_matrix, k_matrix, core_radius_in_um, fibre_radius_in_um');

        best_fitness_id += 1;

    return objective

# The registration has already been performed. Load the results.
if os.path.isfile("outputs/laplacian1.dat"):
    temp = np.loadtxt("outputs/laplacian1.dat");
    sigma_core = temp[0];
    k_core = temp[1];
    sigma_fibre = temp[2];
    k_fibre = temp[3];
    sigma_matrix = temp[4];
    k_matrix = temp[5];
    core_radius = temp[6];
    fibre_radius = temp[7];

# Perform the registration using CMA-ES
else:

    sigma_core = 5.;
    sigma_fibre = 0.75;
    sigma_matrix = 0.6;

    k_core = 1000;
    k_fibre = 1000;
    k_matrix = 1000.0;

    x0 = [
        sigma_core, k_core,
        sigma_fibre, k_fibre,
        sigma_matrix, k_matrix,
        core_radius, fibre_radius
    ];

    bounds = [
        [
            0.005, 0.0,
             0.005, 0.0,
             0.005, 0.0,
             0.95 * core_radius, 0.95 * fibre_radius
        ],
        [
            10.0, 2000,
             2.5, 2000,
             2.5, 2000,
             1.15 * core_radius, 1.15 * fibre_radius
        ]
    ];

    best_fitness = sys.float_info.max;
    best_fitness_id = 0;
    prefix = "laplacian1_";

    opts = cma.CMAOptions()
    opts.set('tolfun', 1e-4);
    opts['tolx'] = 1e-4;
    opts['bounds'] = bounds;
    opts['CMA_stds'] = [0.25, 20.25, 0.25, 20.25, 0.25, 20.25, core_radius * 0.1, fibre_radius * 0.1];

    es = cma.CMAEvolutionStrategy(x0, 0.25, opts);
    es.optimize(fitnessFunctionLaplacian);

    sigma_core = es.result.xbest[0];
    k_core = es.result.xbest[1];
    sigma_fibre = es.result.xbest[2];
    k_fibre = es.result.xbest[3];
    sigma_matrix = es.result.xbest[4];
    k_matrix = es.result.xbest[5];
    core_radius = es.result.xbest[6];
    fibre_radius = es.result.xbest[7];

    np.savetxt("outputs/laplacian1.dat", [sigma_core, k_core, sigma_fibre, k_fibre, sigma_matrix, k_matrix, core_radius, fibre_radius], header='sigma_core, k_core, sigma_fibre, k_fibre, sigma_matrix, k_matrix, core_radius_in_um, fibre_radius_in_um');

(5_w,10)-aCMA-ES (mu_w=3.2,w_1=45%) in dimension 8 (seed=879903, Thu Apr 29 18:44:16 2021)
Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]
    1     10 1.076012272798531e-01 1.0e+00 2.46e-01  6e-02  5e+00 0:44.1
    2     20 1.072723149959538e-01 1.2e+00 2.37e-01  6e-02  5e+00 1:10.3
    3     30 1.074260206848244e-01 1.3e+00 2.33e-01  6e-02  5e+00 1:27.5
    4     40 1.070165335870282e-01 1.5e+00 2.30e-01  5e-02  5e+00 2:02.9
    5     50 1.070399120218351e-01 1.6e+00 2.07e-01  5e-02  4e+00 2:20.4
    6     60 1.071295037961749e-01 1.6e+00 2.00e-01  4e-02  4e+00 2:37.8
    7     70 1.071295465152031e-01 1.6e+00 2.03e-01  5e-02  4e+00 2:55.2
    8     80 1.073673557296478e-01 1.7e+00 2.23e-01  5e-02  4e+00 3:12.5
    9     90 1.071416242594982e-01 1.8e+00 2.40e-01  6e-02  5e+00 3:30.0
   10    100 1.075464775960256e-01 1.9e+00 2.16e-01  5e-02  4e+00 3:47.3
   11    110 1.072859343503126e-01 1.9e+00 2.06e-01  5e-02  4e+00 4:04.8
   12    120 1.069652222743676e-01 1.9e+00 1.96e-01  4e-02  4e+00 4:31.0
   13    130 1.072279359076233e-01 1.8e+00 1.85e-01  4e-02  4e+00 4:48.5
   14    140 1.070567804764945e-01 1.9e+00 1.77e-01  4e-02  3e+00 5:05.9
   15    150 1.070423992230165e-01 2.0e+00 1.83e-01  4e-02  4e+00 5:23.3
   16    160 1.072146211904530e-01 2.1e+00 1.89e-01  4e-02  4e+00 5:40.7
   17    170 1.070371029231238e-01 2.1e+00 2.07e-01  5e-02  5e+00 5:58.1
   18    180 1.071407910789908e-01 2.3e+00 2.28e-01  5e-02  6e+00 6:15.6
   20    200 1.073467092548721e-01 2.4e+00 1.98e-01  4e-02  5e+00 6:50.4
   22    220 1.072555406472326e-01 2.5e+00 1.60e-01  3e-02  4e+00 7:25.3
   24    240 1.070376710759731e-01 2.8e+00 1.41e-01  2e-02  3e+00 8:00.4
   26    260 1.069148812883883e-01 3.1e+00 1.34e-01  2e-02  3e+00 8:44.2
   28    280 1.070435085725500e-01 3.3e+00 1.21e-01  2e-02  3e+00 9:19.0
   29    290 1.068394423876763e-01 3.7e+00 1.15e-01  2e-02  3e+00 9:54.2
   31    310 1.069047350974530e-01 4.2e+00 1.32e-01  2e-02  4e+00 10:29.2
   33    330 1.068767448723697e-01 4.3e+00 1.34e-01  2e-02  4e+00 11:04.2
   34    340 1.067530295837553e-01 4.2e+00 1.47e-01  2e-02  4e+00 11:30.6
   36    360 1.068038409963439e-01 4.6e+00 1.68e-01  2e-02  4e+00 12:06.3
   38    380 1.068231728585128e-01 4.9e+00 1.72e-01  2e-02  5e+00 12:47.1
   40    400 1.066409059532875e-01 5.2e+00 1.77e-01  3e-02  5e+00 13:35.8
   42    420 1.066845977151266e-01 5.1e+00 1.86e-01  3e-02  5e+00 14:19.5
   44    440 1.065885430654243e-01 5.2e+00 1.50e-01  2e-02  4e+00 15:20.8
   46    460 1.065104206693342e-01 5.4e+00 1.49e-01  2e-02  4e+00 16:13.6
   48    480 1.065051023398488e-01 5.9e+00 1.53e-01  2e-02  4e+00 16:57.6
   49    490 1.064991236843169e-01 6.1e+00 1.55e-01  2e-02  4e+00 17:33.2
   50    500 1.064836003924008e-01 6.5e+00 1.48e-01  2e-02  4e+00 18:09.3
   52    520 1.064416378693736e-01 6.9e+00 1.33e-01  1e-02  4e+00 19:03.2
   54    540 1.063350295851035e-01 7.4e+00 1.35e-01  2e-02  4e+00 20:05.6
   56    560 1.061791872810207e-01 7.6e+00 1.62e-01  2e-02  4e+00 20:49.8
   58    580 1.060339350399706e-01 7.7e+00 2.19e-01  2e-02  7e+00 21:34.8
   60    600 1.058692685925371e-01 9.8e+00 2.99e-01  3e-02  1e+01 22:46.2
   62    620 1.053438584504779e-01 1.1e+01 3.97e-01  4e-02  1e+01 23:48.1
   64    640 1.051042075827869e-01 1.2e+01 5.52e-01  6e-02  2e+01 24:33.2
   66    660 1.045866758955183e-01 1.3e+01 6.00e-01  6e-02  2e+01 25:17.1
   69    690 1.048737864699351e-01 1.5e+01 7.95e-01  8e-02  3e+01 26:09.9
   72    720 1.045317618749732e-01 1.7e+01 6.75e-01  7e-02  2e+01 27:11.0
   74    740 1.038606820862940e-01 1.9e+01 6.86e-01  7e-02  3e+01 28:21.3
   77    770 1.035901016947094e-01 2.3e+01 6.69e-01  6e-02  3e+01 29:31.2
   80    800 1.035573198702689e-01 2.6e+01 7.69e-01  6e-02  4e+01 30:32.3
   83    830 1.032688780547188e-01 3.3e+01 7.17e-01  5e-02  3e+01 31:42.4
   86    860 1.032210229467459e-01 3.4e+01 6.88e-01  5e-02  3e+01 32:43.5
   89    890 1.032603936229804e-01 3.7e+01 5.56e-01  4e-02  3e+01 33:35.6
   92    920 1.030591641436899e-01 3.9e+01 6.34e-01  4e-02  3e+01 34:46.1
   95    950 1.030892441052953e-01 4.1e+01 6.06e-01  4e-02  3e+01 35:52.1
   98    980 1.030746888513645e-01 4.3e+01 6.28e-01  4e-02  3e+01 36:59.2
  100   1000 1.029522230627630e-01 4.4e+01 6.31e-01  3e-02  3e+01 37:58.8
  102   1020 1.029472805728889e-01 4.7e+01 5.99e-01  3e-02  3e+01 38:57.1
  104   1040 1.028916785516822e-01 4.8e+01 4.63e-01  2e-02  2e+01 39:58.4
  107   1070 1.029059969994695e-01 4.8e+01 3.46e-01  2e-02  1e+01 41:02.9
  110   1100 1.028649409838554e-01 5.0e+01 2.86e-01  1e-02  1e+01 42:12.2
  112   1120 1.028253978105590e-01 5.0e+01 2.72e-01  1e-02  1e+01 43:13.1
  115   1150 1.028000735858408e-01 5.4e+01 2.27e-01  9e-03  9e+00 44:37.0
  119   1190 1.027726947024870e-01 5.8e+01 2.31e-01  9e-03  8e+00 45:59.8
  123   1230 1.028017244433488e-01 6.9e+01 1.76e-01  7e-03  6e+00 47:15.0
  127   1270 1.027852084406508e-01 7.2e+01 1.93e-01  8e-03  7e+00 48:27.0
  130   1300 1.027651145153484e-01 6.6e+01 1.80e-01  7e-03  6e+00 49:37.3
  134   1340 1.027642440571327e-01 8.0e+01 2.25e-01  9e-03  8e+00 50:56.8
  138   1380 1.027104621467252e-01 9.0e+01 2.72e-01  1e-02  1e+01 52:25.8
  142   1420 1.027871081934873e-01 9.5e+01 2.85e-01  1e-02  1e+01 53:36.4
  146   1460 1.027639314512475e-01 9.9e+01 1.90e-01  7e-03  7e+00 54:48.1
  150   1500 1.027497405832276e-01 1.0e+02 1.77e-01  6e-03  6e+00 55:59.9
  155   1550 1.027466398373111e-01 9.5e+01 1.41e-01  4e-03  5e+00 57:27.6
  158   1580 1.027588067120279e-01 1.0e+02 1.04e-01  3e-03  3e+00 58:20.3

if not os.path.exists("plots/laplacian1_registration.gif"):
    registration_image_set = createAnimation("outputs/laplacian1_simulated_CT_",
                'plots/laplacian1_registration.gif');

Apply the result of the registration

# Load the matrix
setMatrix(matrix_geometry_parameters);

# Load the cores and fibres
setFibres(centroid_set);

gvxr.saveSTLfile("fibre", "outputs/laplacian1_fibre.stl");
gvxr.saveSTLfile("core",  "outputs/laplacian1_core.stl");

print("Core diameter:", round(core_radius * 2), "um");
print("Fibre diameter:", round(fibre_radius * 2), "um");

Core diameter: 16 um
Fibre diameter: 107 um

# Simulate the corresponding CT aquisition
simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram(
    [sigma_core, sigma_fibre, sigma_matrix],
    [k_core, k_fibre, k_matrix],
    ["core", "fibre", "matrix"]
);

# Reconstruct the CT slice
simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);
simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);
normalised_simulated_CT = (simulated_CT - simulated_CT.mean()) / simulated_CT.std();

# Compute the ZNCC
print("ZNCC phase contrast registration:",
      "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, normalised_simulated_CT))));

ZNCC phase contrast registration: 91.66

Optimisation of the phase contrast and the LSF

old_lsf = copy.deepcopy(lsf_kernel);

def fitnessFunctionLaplacianLSF(x):
    global best_fitness;
    global best_fitness_id;
    global prefix;

    global lsf_kernel;

    # sigma_core = x[0];
    k_core = x[0];
    # sigma_fibre = x[2];
    k_fibre = x[1];
    # sigma_matrix = x[4];
    k_matrix = x[2];

    a2 = x[3];
    b2 = x[4];
    c2 = x[5];
    d2 = x[6];
    e2 = x[7];
    f2 = x[8];

    # The response of the detector as the line-spread function (LSF)
    t = np.arange(-20., 21., 1.);
    lsf_kernel=lsf(t*41, a2, b2, c2, d2, e2, f2);
    lsf_kernel/=lsf_kernel.sum();

    # Simulate a sinogram
    simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram(
        [sigma_core, sigma_fibre, sigma_matrix],
        [k_core, k_fibre, k_matrix],
        ["core", "fibre", "matrix"]
    );

    # Compute the objective value
    objective = metrics(reference_sinogram, simulated_sinogram);

    # The block below is not necessary for the registration.
    # It is used to save the data to create animations.
    if best_fitness > objective:
        best_fitness = objective;

        # Save the simulated sinogram
        simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);
        saveMHA("outputs/" + prefix + "simulated_sinogram_" + str(best_fitness_id) + ".mha",
                simulated_sinogram,
                [pixel_spacing_in_mm, angular_step, pixel_spacing_in_mm]);

        # Reconstruct the CT slice
        simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);

        # Save the simulated CT slice
        saveMHA("outputs/" + prefix + "simulated_CT_" + str(best_fitness_id) + ".mha",
                simulated_CT,
                [pixel_spacing_in_mm, pixel_spacing_in_mm, pixel_spacing_in_mm]);

        np.savetxt("outputs/" + prefix + "laplacian_" + str(best_fitness_id) + ".dat", [k_core, k_fibre, k_matrix], header='k_core, k_fibre, k_matrix');
        np.savetxt("outputs/" + prefix + "LSF_" + str(best_fitness_id) + ".dat", [a2, b2, c2, d2, e2, f2], header='a2, b2, c2, d2, e2, f2');


        best_fitness_id += 1;

    return objective;

# The registration has already been performed. Load the results.
if os.path.isfile("outputs/laplacian2.dat") and os.path.isfile(output_directory + "/lsf2.dat"):
    temp = np.loadtxt("outputs/laplacian2.dat");
    k_core = temp[0];
    k_fibre = temp[1];
    k_matrix = temp[2];

    temp = np.loadtxt("outputs/lsf2.dat");
    a2 = temp[0];
    b2 = temp[1];
    c2 = temp[2];
    d2 = temp[3];
    e2 = temp[4];
    f2 = temp[5];

# Perform the registration using CMA-ES
else:

    a2 = 601.873;
    b2 = 54.9359;
    c2 = -3.58452;
    d2 = 0.469614;
    e2 = 6.32561e+09;
    f2 = 1.0;

    x0 = [
        k_core,
        k_fibre,
        k_matrix,
        a2, b2, c2, d2, e2, f2
    ];

    bounds = [
        [
            k_core-500,
            k_fibre-500,
            k_matrix-500,
            a2 - a2 / 4.,
            b2 - b2 / 4.,
            c2 + c2 / 4.,
            d2 - d2 / 4.,
            e2 - e2 / 4.,
            f2 - f2/ 4.
        ],
        [
            k_core+500,
            k_fibre+500,
            k_matrix+500,
            a2 + a2 / 4.,
            b2 + b2 / 4.,
            c2 - c2 / 4.,
            d2 + d2 / 4.,
            e2 + e2 / 4.,
            f2 + f2/ 4.
        ]
    ];

    best_fitness = sys.float_info.max;
    best_fitness_id = 0;
    prefix = "laplacian2_"

    opts = cma.CMAOptions()
    opts.set('tolfun', 1e-4);
    opts['tolx'] = 1e-4;
    opts['bounds'] = bounds;
    #opts['seed'] = 987654321;
    # opts['maxiter'] = 5;
    opts['CMA_stds'] = [1250 * 0.2, 1250 * 0.2, 1250 * 0.2,
        a2 * 0.2, b2 * 0.2, -c2 * 0.2, d2 * 0.2, e2 * 0.2, f2 * 0.2];

    es = cma.CMAEvolutionStrategy(x0, 0.25, opts);
    es.optimize(fitnessFunctionLaplacianLSF);

    k_core = es.result.xbest[0];
    k_fibre = es.result.xbest[1];
    k_matrix = es.result.xbest[2];

    a2 = es.result.xbest[3];
    b2 = es.result.xbest[4];
    c2 = es.result.xbest[5];
    d2 = es.result.xbest[6];
    e2 = es.result.xbest[7];
    f2 = es.result.xbest[8];

    np.savetxt("outputs/laplacian2.dat", [k_core, k_fibre, k_matrix], header='k_core, k_fibre, k_matrix');
    np.savetxt("outputs/lsf2.dat", [a2, b2, c2, d2, e2, f2], header='a2, b2, c2, d2, e2, f2');

(5_w,10)-aCMA-ES (mu_w=3.2,w_1=45%) in dimension 9 (seed=790393, Thu Apr 29 19:43:26 2021)
Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]
    1     10 1.026285346850123e-01 1.0e+00 2.34e-01  2e-02  3e+08 0:43.5
    2     20 1.023785373035976e-01 1.2e+00 2.18e-01  2e-02  3e+08 1:11.2
    3     30 1.024207639063879e-01 1.3e+00 2.04e-01  2e-02  3e+08 1:26.8
    4     40 1.025757418321224e-01 1.3e+00 1.86e-01  2e-02  2e+08 1:43.0
    5     50 1.024522131079754e-01 1.4e+00 1.77e-01  2e-02  2e+08 1:59.7
    6     60 1.024899964278362e-01 1.5e+00 1.66e-01  1e-02  2e+08 2:18.0
    7     70 1.023432758364960e-01 1.6e+00 1.64e-01  1e-02  2e+08 2:42.5
    8     80 1.021204130634706e-01 1.7e+00 1.76e-01  2e-02  2e+08 3:21.9
    9     90 1.022005649875590e-01 2.0e+00 1.86e-01  2e-02  2e+08 3:37.4
   10    100 1.021609788322490e-01 2.1e+00 2.00e-01  2e-02  3e+08 3:52.7
   11    110 1.021738691957654e-01 2.2e+00 2.05e-01  2e-02  3e+08 4:07.9
   12    120 1.020838460415897e-01 2.3e+00 2.15e-01  2e-02  3e+08 4:32.2
   13    130 1.021207578328723e-01 2.5e+00 2.16e-01  2e-02  3e+08 4:47.4
   14    140 1.020464402281892e-01 2.8e+00 2.04e-01  2e-02  3e+08 5:20.7
   15    150 1.019197460463742e-01 2.8e+00 2.21e-01  2e-02  3e+08 5:53.9
   16    160 1.020180723710243e-01 3.2e+00 2.31e-01  2e-02  3e+08 6:09.1
   18    180 1.019194229940256e-01 3.4e+00 2.49e-01  2e-02  3e+08 6:48.6
   20    200 1.019982112309591e-01 3.8e+00 2.43e-01  2e-02  3e+08 7:19.1
   21    210 1.018928413382019e-01 3.9e+00 2.28e-01  2e-02  3e+08 7:43.3
   23    230 1.019226726726832e-01 3.7e+00 2.43e-01  2e-02  3e+08 8:13.8
   25    250 1.018609838117923e-01 3.8e+00 2.19e-01  2e-02  3e+08 8:53.4
   27    270 1.019382790632027e-01 3.9e+00 1.92e-01  2e-02  3e+08 9:23.9
   29    290 1.018736643483591e-01 4.1e+00 1.78e-01  2e-02  2e+08 9:54.4
   31    310 1.018624835093228e-01 4.4e+00 1.68e-01  1e-02  2e+08 10:24.9
   33    330 1.018927047139252e-01 4.6e+00 1.50e-01  1e-02  2e+08 10:55.4
   35    350 1.018929940990494e-01 5.0e+00 1.67e-01  1e-02  2e+08 11:34.8
   37    370 1.019539434957100e-01 5.7e+00 1.60e-01  1e-02  3e+08 12:05.3
   39    390 1.019159219104646e-01 6.3e+00 1.61e-01  1e-02  3e+08 12:35.7
   41    410 1.018616306351648e-01 6.5e+00 1.55e-01  1e-02  2e+08 13:06.2
   43    430 1.019233059927706e-01 7.0e+00 1.35e-01  1e-02  2e+08 13:36.7
   45    450 1.018828026648591e-01 8.5e+00 1.29e-01  1e-02  2e+08 14:07.1
   48    480 1.018787489822530e-01 8.4e+00 1.15e-01  9e-03  2e+08 14:53.0
   51    510 1.019151934316629e-01 8.5e+00 1.12e-01  9e-03  2e+08 15:38.8
   54    540 1.018870847657683e-01 9.6e+00 9.08e-02  7e-03  1e+08 16:24.5
   56    560 1.018567850340121e-01 1.0e+01 8.28e-02  7e-03  1e+08 17:12.7
   59    590 1.019002813156876e-01 9.8e+00 6.53e-02  5e-03  1e+08 17:58.4

if not os.path.exists("plots/laplacian2_registration.gif"):
    registration_image_set = createAnimation("outputs/laplacian2_simulated_CT_",
                'plots/laplacian2_registration.gif');

Apply the result of the registration

# The response of the detector as the line-spread function (LSF)
t = np.arange(-20., 21., 1.);
lsf_kernel=lsf(t*41, a2, b2, c2, d2, e2, f2);
lsf_kernel/=lsf_kernel.sum();
np.savetxt("outputs/LSF_optimised.txt", lsf_kernel);

# Simulate the corresponding CT aquisition
simulated_sinogram, normalised_projections, raw_projections_in_keV = simulateSinogram(
    [sigma_core, sigma_fibre, sigma_matrix],
    [k_core, k_fibre, k_matrix],
    ["core", "fibre", "matrix"]
);

# Reconstruct the CT slice
simulated_sinogram.shape = (simulated_sinogram.size // simulated_sinogram.shape[2], simulated_sinogram.shape[2]);
simulated_CT = iradon(simulated_sinogram.T, theta=theta, circle=True);
normalised_simulated_CT = (simulated_CT - simulated_CT.mean()) / simulated_CT.std();

# Compute the ZNCC
print("ZNCC phase contrast and LSF registration:",
      "{:.2f}".format(100.0 * np.mean(np.multiply(normalised_reference_CT, normalised_simulated_CT))));

ZNCC phase contrast and LSF registration: 92.43

fig=plt.figure();
plt.title("Response of the detector (LSF)");
plt.plot(t, old_lsf, label="Before optimisation");
plt.plot(t, lsf_kernel, label="Before optimisation");
plt.legend();
plt.savefig('plots/LSF_optimised.pdf');
plt.savefig('plots/LSF_optimised.png');

Extract the fibre in the centre of the CT slices

def findFibreInCentreOfCtSlice():
    global centroid_set;
    global reference_CT;
    global cylinder_position_in_centre_of_slice;

    # Find the cylinder in the centre of the image
    cylinder_position_in_centre_of_slice = None;
    best_distance = sys.float_info.max;

    for centre in centroid_set:
        distance = math.pow(centre[0] - reference_CT.shape[1] / 2,2 ) + math.pow(centre[1] - reference_CT.shape[0] / 2, 2);

        if best_distance > distance:
            best_distance = distance;
            cylinder_position_in_centre_of_slice = copy.deepcopy(centre);

    return cylinder_position_in_centre_of_slice;

findFibreInCentreOfCtSlice();

reference_fibre_in_centre = np.array(copy.deepcopy(reference_CT[cylinder_position_in_centre_of_slice[1] - roi_length:cylinder_position_in_centre_of_slice[1] + roi_length, cylinder_position_in_centre_of_slice[0] - roi_length:cylinder_position_in_centre_of_slice[0] + roi_length]));
test_fibre_in_centre      = np.array(copy.deepcopy(simulated_CT[cylinder_position_in_centre_of_slice[1] - roi_length:cylinder_position_in_centre_of_slice[1] + roi_length, cylinder_position_in_centre_of_slice[0] - roi_length:cylinder_position_in_centre_of_slice[0] + roi_length]));

profile_reference = copy.deepcopy(np.diag(reference_fibre_in_centre));
profile_test = copy.deepcopy(np.diag(test_fibre_in_centre));

reference_fibre_in_centre -= reference_fibre_in_centre.mean();
reference_fibre_in_centre /= reference_fibre_in_centre.std();

test_fibre_in_centre -= test_fibre_in_centre.mean();
test_fibre_in_centre /= test_fibre_in_centre.std();

norm = cm.colors.Normalize(vmax=1.25, vmin=-0.5)

fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
plt.tight_layout()
fig.suptitle('Fibre in the centre of the CT slices')

ax1.set_title("Reference image");
imgplot1 = ax1.imshow(reference_fibre_in_centre, cmap="gray",
                     norm=norm);

ax2.set_title("Simulated CT slice after automatic registration");
imgplot2 = ax2.imshow(test_fibre_in_centre,
                     cmap='gray',
                     norm=norm);

comp_equalized = compare_images(reference_fibre_in_centre, test_fibre_in_centre, method='checkerboard');
ax3.set_title("Checkboard comparison between\n" +
              "the reference and simulated images\nZNCC: " +
              "{:.2f}".format(100.0 * np.mean(np.multiply(reference_fibre_in_centre, test_fibre_in_centre))));
imgplot3 = ax3.imshow(comp_equalized,
                     cmap='gray',
                     norm=norm);

plt.savefig('plots/simulated_CT_slice_with_fibres_after_cube_registration.pdf');
plt.savefig('plots/simulated_CT_slice_with_fibres_after_cube_registration.png');

plt.figure()

plt.plot(profile_reference, label="Reference");
plt.plot(profile_test, label="Simulation");
plt.legend();

183.87853883690607
182.36294362794987